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Background: Dust Obscured Star-Formation at z > 1 and Dust in 

the Outer Solar System 
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ABSTRACT 

We provide measurements of the integrated galaxy light at 70, 160, 250, 350 
and 500 /zm using deep far- infrared and submillimeter data from space (Spitzer) 
and balloon platform (BLAST) extragalactic surveys. We use the technique of 
stacking at the positions of 24 fim sources, to supplement the fraction of the inte- 
grated galaxy hght that is directly resolved through direct detections. We demon- 
strate that the integrated galaxy light even through stacking, falls short by factors 
of 2—3 in resolving the extragalactic far-infrared background. We also show that 
previous estimates of the integrated galaxy light (IGL) through stacking, have 
been biased towards high values. This is primarily due to multiple counting of 
the far-infrared/submillimeter flux from 24 //m sources which are clustered within 
the large point spread function of a brighter far-infrared source. Using models 
for the evolution of the luminosity function at z < 1.2 which are constrained by 
observations at 24 fim and 70 /xm, and which are consistent with the results from 
the stacking analysis, we find that galaxies ed, z < 1.2, account for ~ 95 — 55% of 
the extragalactic far-infrared background in the ~ 70 — 500 fim range respectively. 
This places strong upper limits on the fraction of dust obscured star-formation 
at 2; > 1, which are remarkably, below the values derived from the extinction 
corrected ultraviolet luminosities of galaxies. The largest fraction of the total 
40 — 500 fim EBL comes from galaxies between Lir=L(8 — 1000 fim) = 10^^^^^'^ Lq] 
ultraluminous infrared galaxies with Lir > 10^^ Lq contribute only 5%. We use 
the results to make predictions for the nature of galaxies that extragalactic sur- 
veys with Herschel Space Observatory will reveal. Finally, from our constraints 
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on the far- infrared IGL, we provide evidence for the existence of ice mantle dust, 
orbiting the sun at a distance of ~40 AU, which is contributing intensity to both 
the near- and far-infrared background. The presence of this component which 
is a tiny fraction of the zodiacal light emission even at high ecliptic latitudes, 
eliminates any discrepancy between the integrated galaxy light and the diffuse, 
extragalactic background light at all infrared wavelengths. 

Subject headings: cosmology: observations — diffuse radiation — galaxies:evolution 
— galaxies: high-redshift — Kuiper Belt 



1. Introduction 



The extragalactic background light (EBL)E| at far-infrared (FIR) wavelengths, spanning 
a wavelength range of 40 — 1000 fim, is the sum total of radiation from all astrophysical 
processes over cosmic time that is absorbed by dust and re-radiated. It places a strong 
constraint on the fraction of nucleosynthesis and accretion activity that is obscured by dust. 
Decomposition of the EEL into the contribution from individual galaxies as a function of 
redshift, provides a more complete understanding of the rate of build up of stellar mass and 
the growth of black holes over cosmic time and thereby helps identify the epoch at which 
most of the stars and black holes in the present day Universe formed. 



2 gj, 1 



The FIR EBL is thought to have a total intensity between 40—1000 /im of ~25 nWm 
although values as high as ~55nWm~^sr ~^ have been suggested by DIRBE observations 



(IHauser et al.lll998l : iFinkbeiner et al.ll2000l ). Much of the difference arises due to the diffi- 
culty in subtracting the contribution from zodiacal light emission. The zodiacal light con- 
tributes ~30% of the line of sight intensity at 140/um in even a relatively high ecliptic latitude 
(/3 ~ 45°) field like t he Lockman Hole, with a rapidly increasing fractional contribution at 
shorter wavelengths (IHauser et al.lll998l ). Even small errors in the subtraction of the zodiacal 
light propagate through as large errors in the estimate of the FIR EBL, an issue which has 
plagued measurements of the EBL at 60 /im and 100 /im. 

Despite the fact that it is measured with a factor of ~100 less precision than the 
cosmic microwave background, the FIR EBL is arguably the second most energetic ex- 
tragalactic background after the cosmic microwave background which has a total intensity 
of ~1000 nWm~^ sr~^. The intensity of the FIR EBL is comparable, if not in excess of 
the intensity of the optical/near-infrared extragalactic background, whose measurement is 



Also referred to as the Cosmic Infrared Background 



- 3 - 



affected by similar uncerta inties associated with the zodiacal light contribution ( iBernstein 
20071 : iLevenson et al.l 120071 ). Since the optical/near-infrared EBL measures the redshifted 
contribution of unobscured starlight and AGN activity, integrated over cosmic time, the fact 
that the two intensities are comparable, illustrates the importance of dust obscuration in 
quantifying the bolometric luminosity of star-forming galaxies and active galactic nuclei. 

Ever since the direct measurement of the EBL using the Far-Infrared Ab solute Spec 



trorn e ter (FIRAS) and the Diffuse Infrared Background Experiment (DIRBE; iPuget et al. 



19961 : iHauser et al.lll998l : iFixsen et al.lll998l ). various attempts have been made to identify 
the individual sources contributing to the background. The sum of the contribution of 
individual galaxies is termed the integrated galaxy light (IGL). Us i ng ba ckward evolution 
models which fit multiwavelength number counts, IChary fc ElbazI ( 120011 . hereafter CEOl) 
demonstrated that dust obscured infrared luminous galaxies contribute 85% of the 140 fim 
EBL and that 82% of the peak of the FIR EBL arises from galaxies at 2; < 1.5. The 140 fim 
EBL is ~60% of the total FIR EBL in the CEOl models. This implies that about 50% of 
the total FIR EBL arises at 2; < 1.5. 



Elbaz et al.l ( 120021 ) utihzed model spectral energy distributions to estimate the contri- 
bution to the FIR EBL of infrared luminous galaxies individually detected by the ISOCAM 
in a deep 15 /im survey in the Hubble Deep Field- North. They found that the ISOCAM 
LIRGs at 2 < 1 contribute 16=b 5 nWm~^ sr~^ of t he 14 0/im background. More recently, 
by stacking the far-infrared data, iBethermin et al.l (120101 ) showed that galaxies with 24/im 
flux densities >35 /xJy contribute 9.0±1.1 nWm~^ sr~^ of the 160 ^m background. Based on 
photometric redshifts, these galaxies are primarily thought to be infrared luminous galaxies 
at 2; ~ 1 with a tail extending to ~ 2.5 (ICaputi et al.ll2007l ). 



In this context, the recent results from Balloon Large Aperture Submillimeter Telescope 
(BLAST) have been somewhat surpri sing (IDevlin et al.ll2009l : iMarsden et al.ll2009l ). Using a 
technique similar to lDole et al.l (120061 ) . the BLAST team estimated that stacking on sources 
with 24/im flux density greater than 13 /iJy results in an IGL of 8.6±0.6nWm~^ sr~^ at 
250 /im (i.e. the majority of the EBL at this wavelength). Furthermore, the conclusion is 
that 24 /xm detected galaxies contribute ~80-90% of the total FIR EBL between 250 and 
500 //m. The deep 24 fim data which have been used for their an alysis are only sensitive to 
galaxies brighter than 10^° L© at 2 ~ 1 and >1O^^L0 at 2 ~ 2 ( IChary Il2007al ). It is well 
known that there is a whole population of Lyman-break galaxies which are less luminous 
and which appear to have significant dust obscuration based on their UV-slope. Although 
they should be contributing to the EBL, the BLAST analysis seems to suggest that the 
contribution from fainter galaxies is surprisingly small. 



In this paper, we utilize deep, public Spitzer data taken as part of the Far-Infrared 
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Deep Extragalactic Legacy (FIDEL; P.L M. Dickinson) surve}a at 70 and 160 fim and the 
publicly released BLAST maps at 250, 350 and 500 /im to evaluate the consistency of these 
different estimates. We combine direct detections of sources with stacking analysis to provide 
lower limits to the integrated galaxy light. We utilize recent measurements of the far- infrared 
luminosity function out to ^ ~ 1.3, in conjunction with the model spectral energy distribution 
of galaxies based on various spectroscopic and imaging data, to determine the nature of 
sources contributing to the far-infrared background as a function of redshift and luminosity. 
We place constraints on dust obscured star-formation at z > 1 and discuss implications 
for deep Herschel extragalactic surveys. Finally, we propose the existence of a previously 
unknown cloud of high albedo dust grains orbiting the Sun at a distance of ~40 AU which 
contributes significantly to the intensity of the near- and far-infrared background. 

Throughout this paper we assume a standard cosmology with Hq = 71 kms^^ Mpc^^, 
Qm = 0.27 and Qa = 0.73. 



2. The Integrated Galaxy Light from Direct Detections and Stacking 

The integrated galaxy light is the sum of the observed light from individually detected 
galaxies at a particular wavelength. Due to the sensitivity limit of a survey, it is difficult to 
directly detect the dominant population of galaxies contributing at a particular wavelength. 
This is particularly difficult at far-infrared wavelengths, since the sensitivity of far-infrared 
detectors are relatively poor compared to those at shorter wavelengths. Furthermore, due 
to the large point spread function, source confusion dominates the limiting flux of a survey. 

It is possible to assess the total intensity contribution from galaxies below the detection 
limit of a survey using the technique of stacking. This requires knowing the positions of the 
galaxies in the survey area from a deeper, less confused survey. The 24 fim observations of 
the FIDEL flelds with Spitzer provide a particularly useful sample of galaxies. The surveys 
at this wavelength are >150 times deeper in flux density than a deep far-infrared survey. 
By sampling redshifted hot dust and polycyclic aromatic hydrocarbon features, 24 /xm sur- 
veys reveal a sample of galaxies which have signiflcant dust in them and are likely to be 
contributing to the emission at far-infrared wavelengths due to the strong observed corre- 
lation between mid-infrared and far-infrared luminosities of galaxies (CEOl). By stacking 
the far-infrared data at the positions of these 24 /im sources, the background limited noise a 
can be reduced by a/\fN where N is the number of sources that are stacked on. In reality 
however, the flux from the wings of individual bright sources and from the uncertainty in 
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the subtraction of these brighter source results in an additional term, such that the noise in 
the stacked image is: 



conf,i 



N 

with crlgnfi being the variance in the residual map in the vicinity of the i*'^ source after it 
has been fit and subtracted. If the sources are not subtracted from the map, then aconf,i is 
simply the Poisson noise from the flux of the i*'' source which is typically larger than the 
variance in the residual map. Similarly, the extended wings of the point spread function 
results in the stacked flux biased high as we will demonstrate in the next section. Thus, to 
avoid contamination of the stacked flux, stacking should generally be undertaken on maps 
in which the individual sources detected with high significance are fit for and subtracted. 

Since stacking can only reveal average properties of the stacked sample of galaxies, 
selecting a homogeneous sample of objects provide more illustrative results than stacking 
a heterogeneous sample such as those selected at optical/near- infrared wavelengths. For 
stacking at far-infrared wavelengths, 24 //m-selected sources are the preferred choice due 
to their sup erior spatial resol ution, depth and stable image quality compared to longer 



wavelengths (IRieke et al.ll2004[ ). 



Furthermore, extreme care must be taken to avoid counting the flux in the wings of 
individual sources multiple times, an effect which is important at far-infrared wavelengths 
with the large point spread function. Care must also be taken to include the clustering prop- 
erties of the stacked sources in the error analysis which otherwise results in underestimates 
of the uncertainty in the stacked flux. 



2.1. Spitzer Far-Infrared Surveys 



We use three extragalactic fields for our IGL analysis: GOODS-N, ECDFS and Extended 
Groth Strip (EGS); all of which were observed at 24 yum as part of the GOODS Spitzer a.nd/ or 
FIDEL Legacy surveys (P.I. Mark Dickinson). The depths of the Spitzer 24 /xm data varies 
from field-to-field. For consistency in our stacking analyses, we chose a 24 fim depth for 
each field such that the completeness limits were comparable. These flux limits are listed in 
Table [1] The 24 /xm source catalogs are con structed using the pos itions of sources detected 
at shorter wavelengths (3 — 8 fim) as priors (iMagnelli et al.ll2009al ). 
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2.1.1. VOfim 

We utilize a mosaic constructed using the public 70 /im observations of the GOODS-N 
field taken as part o f the FIDEL program in conjunction with deep observations described in 



Frayer et al.l (120061 ) . Since the combined mosaic covers the entire GOODS-N area, there is 
corresponding deep 24 ^m coverage over this mosaic. The total field size with deep coverage is 
174.3 arcmin^. The Spitzer beam size at 70/im is 18" Full Width at Half Maximum (FWHM) 
while the beam size at 24/im is 5.7" FWHM (Table [1]). At this spatial resolution, almost 
all extragala ctic sources are point sources. Using positional priors from the much deeper 



24 /im data, iMagnelli et al.l (j2009al ) have used point source fitting to estimate the 70 /xm 
fiux density of each 24 /xm source. Using extensive simulations, they have characterized the 
reliability and completeness of the FIDEL 70 fim survey. The 80% completeness limit of the 
70 /im survey is 3 mJy although individual sources in unconfused regions can be detected 
down to fainter fiux densities. 

We first subtract all reliably detected sources from the 70 fim mosaic. We do this 
by scahng the 70 fj,m point spread function (PSF) to the fiux density of each source and 
subtracting it off. The 127 reliably detected sources which have 70 fj,m signal to noise ratios 
> 6 and which span a fiux density range of 1.2 — 40.5 m Jy are subtracted from the map. We 
conservatively choose 6a because the noise is correlated in the final mosaics as a result of 
which the measured noise is biased low compared to the true noise. The exact magnitude of 
this effect is difficult to quantify due to the filtering of the individual 70 /im frames before 
coaddition, but we expect it to be a factor of ~2. These individually detected sources 
contribute a total IGL intensity of 1.92±0.2nWm"^ sr~^. We then take the residual map 
with these sources subtracted and stack on the location of the remaining 1945 24 /im sources 
spanning a fiux density range of 20 — 910 /iJy. We cut out a square region of width 84" around 
the position of each source, starting with the one with the highest 24 /im fiux density. We 
then mask out a region of radius 16" around the nominal position of the source in the map 
after including it in the stack. The cutouts are made in decreasing order of 24 /im fiux density. 
The masked region is exactly the size of the aperture used to measure the photometry and 
ensures that the fiux in any particular pixel is not counted multiple times. To stack, we 
take a simple sum of the pixel values from the cutouts although we note that taking a sigma 
clipped sum does not result in any appreciable difference. Due to the mask, the image 
cutouts of faint 24 /tm sources which are in the vicinity of brighter sources will mostly have 
values of zero except for the unmasked pixels. However, since the pixel contribution from 
these faint sources have been included when making the cutout for the brighter source, we 
ensure we are not double counting the flux in any pixel. We obtain a strong detection in the 
resultant stack with an IGL intensity of 1.7±0.7nWm^^ sr^^. After application of a 70/im 
calibration correction factor of 1.15, this implies the total contribution of the 24 /im detected 
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sources to the IGL should be 4.2±0.7nWm ^ sr ^ of which approximately half is from the 
directly detected > 6cr 70 /im sources (Figured]). 

We note that the masking is a critical step in the stacking process. In the absence 
of masking, we would have counted the flux in the wings of brighter (but still undetected) 
sources and attributed this to faint sources. The IGL resulting from stacking the faint 
24 /im sources would then have been 3.2ibl.lnWm~^ sr~^, about a factor of 2 bias in the 
contribution of the faint sources. After adding in the contribution of detected sources 
and the 1.15 calibration correction factor, the resultant IGL would have been found to 
be 5.9±1.1 nWm~^ sr~^. 

The quoted uncertainties in the stacked quantities are derived by performing Monte- 
Carlo simulations whereby the positions of the 24 /xm sources are randomly offset by an 
amount which is between 1 — 3 FWHM. This offset ensures that uncertainties in the stacked 
intensities due to imperfect subtraction of the brighter sources is taken into account. It 
also ensures that the intrinsic large scale clustering of 24 /im sources is preserved in the 
estimation of uncertainties. Using a random set of positions distributed across the entire 
map for the Monte-Carlo will underestimate the uncertainties due to these two effects. The 
70 /im image is then stacked at these offset positions and combined by taking a simple sum. 
The process is repeated a 100 times and the standard deviation of the aperture flux, with 
all appropriate aperture corrections gives the uncertainty in the stacked value. We find that 
these uncertainty values are typically 5 times higher than what one would derive by stacking 
on completely random positions on the source subtracted map. 

The stacked value can change depending on the following effects. As expected, the 
total stacked flux is sensitive to the size of the masking region. If the masking radius and 
correspondingly, the aperture for photometry, is increased to 35", we find that the stacked 
flux decreases such that the total IGL, after adding the contribution of detected sources and 
calibration corrections, is 3.3±0.7nWm~^ sr"^. The fraction of pixels that go into the stack 
is however reduced such that after completeness corrections, we obtain a similar IGL to the 
one estimated with a smaller masking radius (see below). Estimation and subtraction of a 
local background from each cutout before stacking only results in a change of 1% in the total 
IGL. Stacking on the 24 /xm positions without sorting on 24 /xm flux density does not result 
in a significant detection. 

The masking procedure however only includes a fraction of the pixels which correspond 
to faint 24 /im flux sources contributing to the stack (Figure [T]). For example, if all sources 
with 24 /im flux densities above 25 /x Jy are considered, only 37% of the pixels of these sources 
are included in the stack while if the flux limit is increased to 57 /iJy, more than 54% of the 
pixels of the sources are unmasked and can be included. To correct for this, we need to 
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measure the contribution to the stacked flux from sources in bins of 24 fim flux and divide 
the stacked flux and uncertainty of each bin by the fractional number of pixels contributing 
to the stack. We call this a completeness correction^]. The differential contribution to the 
IGL from 24 /xm sources in a flux bin (/i < can be calculated only by stacking on all 
sources brighter than /i in decreasing order of flux, including masking, and then subtracting 
the IGL derived from stacking on all sources brighter than /2. Simply considering sources 
between fi and and stacking on a mosaic with no sources subtracted as has been un- 



dertaken by iDole et al.l (120061 ) overestimates the flux although the exact overcounting factor 
depends on the relatively uncertain clustering properties of the faint sources around bright 
sources. We find that stacking on the 838 sources with S24 > 74 /iJy, after application of the 
completeness corrections shown in Figure [1] and appropriate calibration corrections account 
for 4.9±0.5nWm~^ sr~^. The 1107 sources between 20 — 74 yuJy only have one fifth of their 
pixels included in the stack due to the masking. The differential contribution to the stacked 
flux in this bin needs to be corrected upward by a factor of 5 yielding 2.1±0.9nWm~^ sr~^. 
Thus, our best estimate for the IGL from the data, including detections, masking, com- 
pleteness corrections and calibration corrections is 7.0±1.0 nWm~^ sr~^. We note that the 
24 /im catalog itself is 85% complete in the faintest flux bins but we neglect application of 
this second order completeness correction which would add 0.4nWm~^ sr^^ to the best IGL 
estimate derived above. 

An alternate technique that we adopt to measure any bias in our stacking technique is by 
simulating the data. We input point sources into a simulated 70 /im map at the locations of 
the 24 /im sources that we stack on. The 70 /im/24 /zm flux ratios for these simulated sources 
are derived from th e CEOl model t emplat es and the evolution of the infrared luminosity 



function derived by iMagnelli et al.l (j2009al ). We also make sure that this map does not 
contain any 6a sources so we are trying to simulate our residual map which we use for 
stacking. We then repeat the stacking procedure on this simulated image using no masking 
and compared the stacked flux with the total input flux. We find that the stacked flux exceeds 
the input flux by a factor of 1.2. Thus, using the simulations, we estimate that the IGL in 
the maps is probably the sum of 1.92 (detected sources) and 3.2/1.2, which after calibration 
correction factors implies 5.3±1.1 nWm~^ sr~^although the exact value may be susceptible 
to the redshift distribution of the sources and the distribution of their 70/im/24/im flux 
ratios. 

Our analysis therefore suggests that the true 70 /xm IGL from 24 /im galaxies brighter 



■^There is a second completeness correction associated with the incompleteness of the 24 /xm catalog at 
the faintest flux densities. Since we dont know the locations of those sources, we cannot include them into 
the simulations 
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than 25 /iJy is in the range 5.3 — 7.0 nW m ^ sr ^, about 30% of which is due to galaxies that 
are individually detected in the 70 /im map. 



A comparison with previous estimates of the 70 fim IGL reveals that iDole et al.l (|2006[ ) 
estimated an IGL of 5.93±1.02 nWm~^ sr^^ for source with S24 > 60;uJy. By adopting a 
total EBL estimate of 6.4nWm~^ sr~^, they claim that the bulk of the IGL is resolved. We 
next a ttempt to r e produ ce these measurements using an identical set of steps as described 
in the iDole et al.l ( 120061 ) paper. We divide the 24 /im source list into flux bins and stack 
on the sources in each bin separately, with no masking. The stacking is perforr ned on the 
maps which still have the individually detected sources in them as was done by iDole et al. 
( 120061 ). We then add the contribution of the stacked measurement in each 24 //m flux bin 



and obtain a 70 /im IGL value of 7.1±1.4 nWm~^ sr~^ for S24 > 60/xJy (Figured]). This is 
a significant overestimate since the stacking on the sources in the fainter flux bins includes 
sources which fall within the PSF of brighter sources. As a result, the flux in those pixels 
is counted multiple times. If we mask the sources as described earher, we derive a value of 



3.3=bl .O nWm ^ sr ^for sources with S24 > 60/iJy which is significantly below the lDole et al. 



( I2OO6I ) estimates. However, after completeness for the fraction of pixels missed and calibra- 



tion corrections, this number rises to 5.4±0.5n Wm ^sr ^, as disc ussed above. Although the 



difference in the derived value compared to the iDole et al.l (120061 ) estimate is not large, it is 
clear that a systematic bias is introduced depending on how the stacking is done. As we will 
demonstrate in the subsequent sections, this systematic bias gets larger with a larger point 
spread function. Thus, we believe that each of the steps: source subtraction, stacking on 
the residual maps with masking and simulations of the bias need to be performed to obtain 
a reliable stacked result. 



Prayer et al.l (120061 ) estimated the IGL by integrating over the 70 /xm source counts 
distribution in GOODS-N, without stacking. They found that the contribution of directly 
resolved sources with S70 > 1.2mJy is 4.3±0.7nWm~^ sr^^. This is consistent with our 
values derived frorn stacking since it is simply an integral over a source counts distribution. 
Prayer et al.l ( 120061 ) derive a total EBL of 7.4±1.9nWm~^ sr~^by adopting a power-law for 
the differential counts and extrapolating the source counts to zero flux. 



The luminosity function of iMagnelli et al.l ( l2009a( ). if evolved between < 2; < 1.3 



predicts a 70 /im IGL from z < 1.3 galaxies of 9.2nWm ^sr^. This is using the CEOl 



far-infrared templates which has been shown by iMagnelli et al.l ( l2009al ) to be a good fit to 
the 24/im/('l-|-z) and 70//m/(l+z) lu minosities of galaxies. Thus, we conclude that the total 
EBL estimate of iPrayer et al.l (120061 ) is a significant underestimate. As we will discuss later, 
evolutionary models predict that the fraction of the 70 /im EBL which can be attributed to 
galaxies with S24 > 25/iJy is ~80% with a total EBL of 9.5±0.3 nWm~^ sr"^ Since our mea- 
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sured IGL from stacking is between 5.3 — 7.0 nWm ^sr ^, the deepest current observations 
at 70 /im have resolved about ~56-75% of the EBL. 



2.1.2. 160^m 

At 160 ^m, we utihze the pubhcly released first epoch FIDEL data in the Extended 
Groth Strip. We analyze only the deepest 914arcmin^ which has a total exposure time of 
>500s at 160 yum. The corresponding coverage in the 24 /im map is >7200s. Details of the 
160 /im mosaic are in Table [H Starting with the 3.6 /im catalog, we fit for t he fluxes of the 



24 /im sources using the same technique as above and which is described in iMagnelli et al. 



(j2009al ). We then apply a flux density cut of 100 /zJy to the 24 /xm catalog. The coordinates 
of these brighter 24 /xm sources are provided as input to the point source fitting in the 160 /im 
image. The uncertainty in the 160 /im photometry of each source is estimated as the sum 
of the squares of the residuals in the 160 /im map after all sources are fit and subtracted. 
We apply a signal to noise cut of >5(t and 20mJy in the 160 /xm catalog for reliability. 
We find that these reliably detected 160 /im sources contribute 1.4±0.1 nWm~^ sr^^. After 
subtracting out these sources, we stack at the positions of the remaining 24 /xm sources 
following the guidelines described earlier for the 70 /xm stacking. The public release of the 
24 /xm EGS observations are relatively shallow compared t o the GOODS-N d ata described 



earlier, but still deeper than the S'pzteer observations used by lDole et al.l ( l2006l ). We estimate 
the 24 /im data in EGS is complete at 55 /i Jy depth. 

Stacking on sources detected with high confidence in the 24 /xm map, including mask- 
ing, results in a stacked 160 /xm intensity of 1.9±0.3 nWm~^ sr~^implying a total 160 /im 
IGL intensity of 3.3±0.4 nWm~^sr~^. We note that the masking radius used (Tabled]) 
is large compared to the FWHM of the PSF due to the fact that the PSF in the broad 
160 /im filter is more sensitive to the intrinsic spectrum of the source in the sense that 
red sources have broader PSFs than blue sources. To avoid missing flux, we use a larger 
aperture and correspondingly larger masking radius. If we had chosen not to mask the 
sources after including them in the stack, we would have derived a total stacked intensity of 
11.8±1.6nWm~^ sr~^and a total IGL of 13.2±1.7nWm~^ sr~^after including the contribu- 
tion of detected sources (Figured]). 

As we did for the 70 /im data, we need to apply completeness corrections for the pixels of 
the stacked sources that have been masked due to their presence near a brighter 24 /im source. 
The fractional incompleteness in the 160 /im stacking at the same 24 /im flux density is higher 
than at 70 /im since the point spread function is larger (Figured])- We flnd that the fraction 
of pixels included in the stack is 39% if all sources brighter than S24 > 212/iJy are stacked 
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together. This number increases to 52% at 275 /iJy. Thus, after completeness corrections 
to the stacking analysis from sources with S24 > 275^Jy and including the contribution of 
individually detected sources, we can account for 2.9±0.3 nWm~^ sr~^of the 160 fim IGL. We 
cannot account for the contributions of fainter 24 /im sources between 20 — 275 ^ Jy since only 
1.7% of the pixels from these sources are included in the stack. The resultant completeness 
corrections are large and unreliable although the stack from this small fraction of sources 
that can be stacked on does result in a significant detection. 

We can however estimate a lower limit for their contribution by stacking down to the 
maximum flux limit /i such that the contribution of sources between /i and 275 /iJy re- 
sults in a significant stacked detection. We find that in the 24 /xm flux density range of 
125—275 /iJy, we can get a measurement of the stacked flux although with a completeness 
correction factor of ~14. The resultant contribution to the IGL from sources in this bin, 
after completeness corrections, is 4.2±0.9 nWm~^ sr^^. Thus, based on the stacking analysis 
as much as 7.1±1.0 nWm~^ sr^^of the 160 ;um IGL may be resolved. 

As we did for the 70 /im data, we attempt to simulate the 160 /xm map to assess how large 
the bias is if we stack on the residual map with no masking. We input point sources using the 
160 /xm PSF at the locations of the 24 /zm sources that we stack on. The 160 /im/24 fim flux 
ratios for these simulated sources are derived from t he CEOl model templa tes and the evolu- 



tion of the infrared luminosity function derived by iMagnelli et al.l ( l2009al ). We then repeat 
the stacking procedure on this simulated image using no masking and compared the stacked 
flux with the total input flux. We find that the stacked flux exceeds the input flux by a factor 
of 1.4. Thus, using the simulations, we estimate that the IGL in the maps is probably the 
sum of 1.4 (detected sources) and 11.8/1.4, which implies 9.8±2.0 nWm~^ sr~^although the 
exact value may be susceptible to the redshift distribution of the sources and the distribution 
of their 160 /im/24 /im flux ratios. 

We conclude that the measured 160 /xm IGL derived from the deepest Spitzer maps is 
3.3nWm^^sr^^ although simulations and estimates of com pleteness at the fa int flux bins 
yield a value between 7.1—9.8 nWm~^sr^^. In comparison, iDole et al.l (|2006[ ) estimate an 



IGL of 10.7±2.3 by stacking on sources with S9 4 > 60 fiJj from their shallower imaging data 



We note that if we follow the iDole et al.l (120061 ) stacking prescription in this field, we obtain 



a contribution of 15.0±2.0 nWm~^ sr~^from S24 > 60 /iJy sources. The overestimation in 
the stacking derived IGL at 160 /im due to the overcounting of flux in individual pixels is 
illustrated in the plots shown in Figure [H 

We also note that due to the large point spread function, it is not possible to reliably 
measure the contribution to the 160 fim EBL from sources below S24 < 275 /iJy. At this flux 
level, less than 50% of the pixels from faint 24 /im sources contribute to the stacked flux 
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and the differential contribution from sources between 120 — 275 /iJy cannot be measured to 
better than 5a. Higher resolution data with Herschel will enable the contribution of fainter 
sources to be measured as described later. 



The total EBL at 160 fim derived by lFixsen et al.l (119981 ) using FIRAS data is 13.6±0.4 
nWm~^ sr~^ while our models described later yield a value of 14.1±0.7nWm~^ sr~^. Thus, 
the current observations have not resolved the 160 /xm background. 



2.2. BLAST Submillimeter Survey 

BLAST is a 1.8m submillimeter telescope on a balloon, with detectors at three wave- 
lengths: 250, 350 and 500 ^m (Pascale et al. 2008). BLAST completed a 11 day flight 
from Antarctica in December 2006 during which it took very deep images of the Extended 
Chandra Deep Field South (Devlin et al. 2009). 

The BLAST data were released to the public in April 2OOC0. This release includes 
several versions of the signal and noise maps at all three wavelengths. The "decon.fits" files 
(hereafter referred to as the decon maps) have had a whitening filter applied to suppress the 
large scale frequencies - this affects the PSF and causes it to have negative structure and zero 
mean. The "smooth. fits" files have been convolved with the PSF in order to estimate the 
point source fluxes. Convolving with the PSF is equivalent to a minimum fitting of the 
PSF with the data (see Serjeant et al. 2003). This method is generally used to extract point 
sources in submm surveys as it increases the signal to noise ratio (SNR) and aids detection 
of faint (> 3cr) sources. However, this technique works under the assumption that a single 
PSF is a good fit to the data and may not be valid when there are blended and confused 
sources (it is only valid if there is << 1 source per beam, Serjeant et al. 2003). 

Point source catalogs down to 3a from these BLAST maps are given in Devlin et 
al. (2009). Information about the PSF and pixel scale for the BLAST maps is given in 
Table [D 



2.2.1. 250, 350 and 500 fim 

Unlike in the infrared, where stacking is done by excising cutouts of the image, coadding 
the cutouts and doing aperture photometry, in the submillimeter, stacking is usually done 



"'http://blastexperiment.info/results.php 
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on the smoothed (beam-convolved) maps where single pixel values, whi ch contain the tota l 



flux of each source, are averaged to produce the stacked average (e.g. iDevlin et al.ll2009l ). 
This is primarily because of the complex beam shape which can have significant sidelobes. 
The sidelobes can contaminate the aperture photometry of other sources, resulting in a 
systematic uncertainty to the aperture flux which is dependent on the flux distribution and 
separation of the sources (See Appendix for details). 

Although stacking on the BLAST maps has been pre sented in three separate papers 
( jPevlin et al.ll2009t iMarsden et al.ll2009l : iPascale et al.ll2009l ). there are certain improvements 
which can be made to the analysis presented therein. In Marsden et al. (2009) it is argued 
that stacking essentially boils down to taking the covariance of the map with the catalog. 
This simplification is true under the assumption that galaxies are not clustered. We know 
that IR luminous galaxies are clustered more strongly than optically selected field galaxies 
( iGilli et al.l 120071 ). Although the clustering of IR luminous galaxies has been measured in 
redshift bins, it is the clustering strength of faint IR sources in the vicinity of brigher IR 
sources that is difficult to measure. We find that the surface density of 24 /im detected 
sources within a radius R (Tabled]) of a 500 /im source brighter than 20 mJy is twice as high 
as the average surface density of 24 /im detected sources across the entire field (Figure [21] 
in the Appendix). Results from stacking are s ystematically biased high by this clu s tering . 
As a result, the sta cking method employed by iDevlin et al.l (|2009[ ): iMarsden et al.l ( 120091 ): 
Pascale et al.l (120091 ) will result in values biased high since the flux from the many clustered 
sources within each BLAST beam will essentially be counted more than once. Serjeant et 
al. (2008) estimated this over-counting of the stacked flux assuming a clustering strength 
and showed that it was significant even at the much smaller beam of SCUBA on the JCMT. 
There are severals steps that can be taken to test and account for this effect of clustering on 
the stacked signal which we apply to our own BLAST stacking analysis below. 

We perform our BLAST stacking analysis on both the decon maps where we stack 
cutouts and the smoothed maps where we stack single pixel values. We find consistent 
results between these two approachs. Values quoted in this paper are all for stacking image 
cutouts on the decon maps. 

As described earlier for the 70 /im and 160 ii m stacking analysis, we start off with the 
list of 24 /im sources cataloged for the ECDFS by iMagnelli et al.l (l2009al ). Building on the 
techniques presented for the Spitzer far-infrared data, we take three additional steps which 
many previous stacking analysis in the submillimeter have not accounted for: 



We fit the bright detected submillimeter sources and subtract their flux from the map. 
This residual map will then be used for the stacking analysis. The bright point sources 
leave a substantial footprint on the map which causes stacked fluxes to be overestimated 
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if they are left in the map. The contribution to the background from the detections 
will be counted separately after doing the stacking analysis. 

• We mask regions of the map as they are included in the stack. Again the large beam 
means that there are many 24 fim sources per beam. If we don't exclude regions of the 
map that are already in the stack then flux in those regions will be counted more than 
once. 

• We perform simulations of confused maps with different source count distributions to 
assess the presence of any bias in the stacked intensity values. 

We calculate the uncertainties on all stacked quantities by doing Monte Carlo simu- 
lations stacking random positions in the map. We performed 10000 trials of stacking 
galaxies (where is the number of galaxies in the catalog). The distribution of stacked 
values from the 10000 trials is confirmed to be Gaussian centered around zero and we take 
the standard deviation of the distribution as our uncertainty in the stacked flux. 

Subtracting the suhmillimeter detections 

In order to test the effects of bright sources on the stacked signals, we iteratively fit the 
sources out of the map. Starting with the brightest 250 /um source, we fit the actual PSF 
(including positive and negative structure) to the decon image at the source position and 
remove the signal from the decon map. We continue to the next brightest source and so on. 
For the positions of the detected sources we use the source lists from Devlin et al. (2009) down 
to 4cr. After removing all the 250 //m emission from detected sources above a given SNR 
threshold, we are left with a residual map (see Figure E]). The distribution of pixel values in 
this residual map compared to the input map is shown in Figure [20] in the Appendix. We 
stack this map on the 24 /xm galaxies. 

We calculate the total flux removed from the maps from the detected sources. Table 
H] lists the contribution to the background from sources detected above 6, 5 and 4(T in the 
central region of the BLAST maps that overlaps the deepest 24 /im coverage. 

We note that there is no bias associated with estimating the contribution from the 
detected sources separately from the stacked sources based on the covariance argument in 
Marsden et al. (2009). The only complication associated with the subtraction of sources 
from the maps might be flux boosting but since source subtraction is limited to the highest 
SNR sources and we are fitting the full PSF (and not simply scaling to the peak flux), flux 
boosting is expected to be much less of an issue. 
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The results of stacking on the residual maps are shown in Table [3] and Figures H] (bar 
plots). Leaving the detected sources in the maps prior to stacking, results in a noticeable 
over estimate of the stacked flux as was shown for the 70 and 160 //m data. By simply 
subtracting the submillimeter detections prior to stacking we found that the IGL resolved 
by 24 /zm-selected galaxies is < 50% the total EBL at these submillimeter wavelengths. 
Stacking on the residual maps will still result in an overestimate since the flux surrounding 
sources < 4cr can still be double counted. We can't subtract these sources individually but 
we can mask them so their flux is not counted multiple times. 



Masking while stacking 

We stack image cutouts from the BLAST maps centered on each 24 fim source which is not 
matched to a detected submillimeter source. Once a cutout is made from the map and added 
to the stack, these pixels are masked out in the original map so that the same flux cannot 
be counted more than once. We step through the 24 /xm catalog in order of decreasing 24 
lim flux since the brighter 24/im sources are expected to have higher 250/im flux on average. 
We ensure that the mean of the map is zero prior to beginning our stacking analysis which 
means that we do not need to do any additional sky-subtraction from our stack. 

We tried a range of cutout sizes and settled on the values given in Table [1] since those 
require a small aperture correction and doesn't mask out to o much additional flu x. Assuming 



a perfect Gaussian beam with FWHM of 36" at 250 /im (IMarsden et al.l l2009l ) this cutout 
size should include all flux within ±3 times the standard deviation of the Gaussian (in reality 
the beam is not a perfect Gaussian). 

We measure the stacked flux from the cutout by fitting the actual PSF. This is better 
than fitting a simple Gaussian or doing aperture photometry since the negative sidelobes 
in the PSF is significant. We tested this on cutouts of detected 250 /im sources and found 
that it gave better results than aperture photometry or fitting a Gaussian (see Figure [T9|) . 
From this test, we derive a negligible aperture correction of 1.02 to account for the emission 
outside the cutout 

When masking while doing the stack, we are only able to stack a small fraction of the 
24 /im sources. Since we are stacking in order of decreasing 24/im fiux we can find the fiux 
at which there is no regions of the map left to stack. The 24/im fiux below which we begin 
to be incomplete is 600/tJy (see Figure [3]). We note that this fiux limit is comparable to that 
obtained at 160 /im (see Fig. [T]) which is consistent since the beams are almost the same size. 
We also note that this is the fiux limit where there are ~6 250/im BLAST beams/source 
which is well past the standard confusion limit of 40 beams/source and even past the 7 
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beams / source limit where pri or based source extraction techniques start to give inadequate 
results dMagnelh et"aDl2009al ). 

As we did for the 70 and 160 /im analysis, we explored applying completeness correc- 
tions to the bins below 600 /iJy. The stacked fluxes for the individual bins are however not 
signiflcant and so we can't reliably apply a correction to account for the sources we have 
had to mask. In a very wide area survey, where there are 100s of objects in each bin below 
600 /iJy then we could imagine applying a completeness correction to the masked stacking 
values to resolve more of the background. This will be possibly with Herschel where the 
beam is half the size so many more sources can be stacked with masking. 

Our results from stacking with masking are shown in Table [3] and Figures HI Since we 
are unable to apply completeness corrections, these values are lower limits to the true IGL 
from the 24 /im galaxies; the true value will lie somewhere in between the masking values 
and the no masking values. In addition, we expect the stacking with masking values will 
be biased low since we have only masked the positive emission in the zero-mean BLAST 
PSF and therefore are stacking on a map with a slight negative mean. From our masking 
analysis, we conclude that with the spatial resolution of the current data, the lower limit to 
the IGL from 24 /im selected galaxies between 250 — 500 /im is ~10-15% of the total EBL 
shown in Table 5. 



Simulations 

We can further test how the stacking analysis might be biased by performing simulations 
where we know the IGL we put into the map (IGLj„) and then we perform the same stacking 
procedure to see what we measure (IGLout)- The ratio of these two quantities is then our 
over-counting or bias factor. This can give us an estimate of the true IGL coming from the 
24 //m-selected galaxies. 

We performed simulations of the BLAST 250 /im maps by doing the following (the same 
steps were repeated for the 350 and 500 /im maps as well): 

1) Assign a 250/im flux to each 24 fim source using the full range of conversions from 
the CEOl galaxy templates. These sources are inserted into the fake map using the actual 
BLAST PSF at the real 24 fim positions so as to include the inherent clustering of this 
population in our simulation. 

2) Add noise to the fake map according to the actual BLAST noise maps. 

3) Subtracted any 4cr detections from the maps to create a fake residual map for stacking. 
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4) Perform stacking at positions of all input galaxies to get lGLout,nomask- 

5) The ratio of lGLin/lGLout,nomask tells us how much we are over-predicting the stacked 
flux when we stack on the residual maps without masking. This bias factor is listed in Table 
il 

6) We can also do the stacking with the masking to get IGLout.mask and the ratio 
lGLin/lGLout,mask tcll US how much we are under-predicting the stacking flux when we mask. 

7) Iteratively determine the IGLj„ which reproduces both our observed lGLout,nomask 
and our observed IGLiM^mask- This value is referred to as IGLcorr- 

The value of IGLcorr is an estimate of the actual background produced from the 24 /im 
galaxies (assuming all 250 fim detections have 24 fim counterparts) and is plotted as the 
blue circle in Figures H] and given in Table HI 

We note th at our values aresignificantly below the FIRAS determined EBL values 



as measured by iFixsen et al.l ( 1l998l ) and therefore conclude that ~45-55% of the submil- 



limeter background curr ently remains u r iresol ved. This conclusion is different from that of 



Devlin et al.l ( 120091 ) and iMarsden et al.l (120091 ) who do not do any masking or subtract the 



submillimeter detections prior to stacking. We have demonstrated that these two steps are 
crucial to ensure that the results are not biased due to overcounting the flux from multiple 
sources within the large beam because of the clustering of IR luminous galaxies. Aside from 
differing stacking techniques we also show later in this paper that we do not expect 24 /im- 
selected galaxies to resolve all the EBL since there are significant amounts of star formation 
coming frorn less lu minous galaxies. We also note that the "completeness" corrections that 



Devlin et al.l ( l2009l ) apply to the FIDEL 24 ^um catalogs are small but incorrect. The com- 
pleteness corrections they apply are based on much shallower data than the FIDEL data 
and the prior based cataloging technique used to generate the FIDEL 24 /im catalogs are 
less affected by confusion, resulting in higher completeness. 

This bias factor listed in Table S] provide a rough guide as to how much the stacked 
signal will overestimate the true value, if one does not correct for the contribution of multiple 
sources per beam due to the clustering of IR luminous galaxies. If all things were equal then 
this bias factor would only depend on the FWHM of the beam and the density of sources 
that are being stacking. However, in practice the bias factor will also vary depending on the 
wavelength since the dominant sources of the background change as a function of wavelength, 
on subtle differences in the PSF compared to a perfect Gaussian and on the clustering 
strength of the sources being stacked. For these reasons, these bias factors cannot be applied 
blindly to other far-infrared and submillimeter stacking analyses. 
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We repeated the above simulations with unclustered, random distribution of the 24 fxm 
sources and found a bias factor of 1.00±0.06. This shows that the dominant effect in estimat- 
ing the stacked fluxes of sources below the detection limit indeed arises from the clustering 
properties of the sources. If they were unclustered, there would not be a bias. 

We also performed the above simulations without step 3 to test stacking biases without 
first subtracting the detections (as Devlin et al. 2009 and Marsden et al. 2009 have done). 
As expected, we find even higher bias factors than those listed in Table 4, with a value of 1.6 
at 250 fj,m instead of 1.16. Factoring in this bias factor at 250 /im, we get IGLcorr = 4.6 ±0.4 
nWm~^ sr~^ which is consistent with our best estimate of IGLcorr = 5.2 ± 0.5 nWm~^ sr~^ 
from the 24 /um galaxies. 



2.3. Cosmic Variance 



Although the submillimeter stacking can only be undertaken with the deep BLAST 
observations taken in the ECDF-S, we can compare the stacking values at 70 fim between 
the GOODS-N and ECDF-S observations to assess if cosmic variance plays a significant role 
in different measures of the IGL. We choose only the central 777arcmin^ of ECDFS which 
has 70 fim coverage of at least half its peak exposure time. Using a procedure identical 
to that described earlier in Section 2.1.1, we find that the 179 individually detected > 6a 
sources in the 70 ;um map contribute 1.72 nWm~^ sr^^. Stacking on the remaining 6153 
24 yum > 5cr sources in this region, down to a flux limit of 25 /i Jy, yields a stacked intensity 
of 1.2±0.4nWm~^ sr~^. After the factor of 1.15 calibration correction, we find a total IGL 
intensity in ECDF-S of 3.1±0.4nWm~^ sr~^. If we had done the stacking without masking, 
we would have obtained a value of 2.4nWm~^sr~^ for the stacked intensity from the indi- 
vidually undetected sources, similar to the factor of 2 bias that we obtained in GOODS-N. 
Although there seems to be cosmic variance between the fields which results in a 26% lower 
IGL in ECDF-S than in GOODS-N, we conclude that our method gives consistent values for 
the IGL in the two fields. 



3. The Evolution of the Infrared Luminosity Function 



The evolution of the far-infrare d luminosity function of star-forming ga l axies has been 
meas ured very accurately at z < 1.3 ( Magnelli et al.|[2009a : Caputi et al.ll2007 : Le Floc'h et al, 
20051 ) . For each redshift, one can adopt a spectral energy distribution for galaxies as a func- 
tion of far-infrared luminosity and calculate the observed intensity at each wavelength from 
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galaxies in each redshift bin. The sum of these observed intensities, integrated over all red- 
shifts and luminosity, is the total EBL. By integrating between < z < 1.3, it is trivial to 
estimate what fraction of the EBL is directly attributable to galaxies at these redshifts and 
if that is consistent with the fraction of EBL that is resolved through stacking. We neglect 
the contribution of active galactic nucl ei (AGN) since it has been demonstrated that their 
net contribution to the EBL is < 10% (ITreister fc Urryl 120061 ). 



The local IR LF is commonly represented as a double power-law (jSanders et al.ll2003l ). 
This LF has a shape of: 
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MagneUi et al.l fl2009af ) use the deep 24 and 70 fim imaging data in the GOODS/FIDEL 



fields, in conjunction with spectroscopic and photometric redshifts to estimate the redshift 
evolution of the 24/xm/(H-z) and 70/im/(l-|-z) luminosity function. By using a combination 
of direct detections and stacking at 70;um the luminosity function can be directly measured 
to fainter luminosities than before, thereby placing the first constr a ints on the faint end slope 
of the infrared luminosity function at z ~ 1 — 2. iMagnelli et al.l ( l2009al ) fit model spectral 
energy distribution to the 24yum and 70/im luminosities of galaxies to estimate the total IR 
luminosity function. They have also shown that the CEOl templates provide an excellent fit 
to the observed 24 and 70 /im luminosities of galaxies at these redshifts. They find that the 
IR luminosity function shows rapid lumin osity evolution as (l-t- z)^'^"*^"'^ between < z < 1 
but appears to flatten at higher redshifts (iMagnelli et al.ll2009bl ). 



We attempt to constrain the evolution of the infrared luminosity function at 2; > 1.3 
using the measured EBL as the primary constraint. In order to calculate the resultant EBL 
arising for different evolutionary histories for the infrared luminosity function, we need to 
adopt a far- infrared spectral energy distribution. CEOl provided a set of model spectral 
energy distribution for infrared luminous galaxies as a function of mid-infrared luminosity. 
Since the mid-infrared luminosity is correlated with the FIR luminosity, we can associate 
these templates with each luminosity bin in the LF to make extrapolations to the other 
wavelengths. 
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3.1. 



Constraints on Dust Enshrouded Star-Formation at 2; > 1: Modified CEOl 

Templates 



We adopt the parametric evolution of iMagnelli et al.l (l2009bl ) whereby (po evolves with 
redshift as (1 + z)"° while Lq evolves with redshift as (1 + z)"-^. aL=3.5 between < z < 1 
and ao = —0.9 between < z < 1. At higher redshifts, we leave the evolution as a free 
parameter and assess the range of evolutionary parameters that are consistent with the FIR 
EBL measured by FIRAS. 

We adopt the luminosity dependent CEOl SEDs at z < 1.3. At 2; > 1.3 we apply the 
(1 + zY'^ luminosity evolution observed between < 2; < 1 to the CEOl SEDs. For example, 
aX z > 1.3, we assume that a Lm = 10^^ galaxy has the SED of a 10^^ x 2.15~^'^L0 in 
the local Universe. The reason we use 2.15 as the multiplicative factor is because the highest 
redshift bin in which the evolution is measured is 1 < z < 1.3 and the mid point of the bin 
corresponds to a 1 + 2 = 2.15. So, we scale the lower luminosity galaxy up by the ratio of 
the luminosities to represent a high redshift galaxy of a higher luminosity. We refer to these 
as the modified CEOl templates. 

We find that at z > 1, if we ke e p Q!r,= 1.3, which is the tentative evolution measured 
at 1 < 2 < 2.3 by iMagnelh et aD J2009bl ). we find that -3.2 < an < -2.3 to avoid 
overproducing the FIRAS measured EBL. This is shown in Figure O We also find that the 
luminosity function evolution measured at 2; < 1.3 accounts for ~ 80 — 90% of the total IGL 
between 60 — 160 /um. This number drops to ~67% at 250 /xm, 55% at 350 //m and 44% at 
500 /um. The statistical uncertainty on these percentages are large since the absolute level of 
the background at 500 /xm is not known to better than 30% (Table 5). 



We note that there exists a discrepancy betwee n the FIRAS EBL (IFixsen et al.lll998l ) 
and the DIRBE measure EBL (IHauser et al.l Il998l ) at A << 200 /im. This difference is 
inconsequential to the evolution at z > 1 as can be seen by the range of EBL values spanned 
by the different star-formation rate densities in Figure O The evolution at 2; > 1 primarily 
affects the EBL values at wavelengths of Xpeak x (1 + -2) where Xpeak is the peak in the rest- 
frame far-infrared SED (z^L^) of galaxies. Xpeak is 60 — 70 um for starbursts l ike M82 and 
even for low-metallicity blue compact dwarf galaxies lEngelbracht et al.l (120081 ). Thus, any 
effect on the EBL from evolution between 1 < 2; < 6 would be at A ~ 120 — 500 /im as has 
been discussed previously in CEOl. 



In an alternate scenario, we adopt the same luminosity function evolution as lMagnelli et al, 
(j2009bl ) at 2; < 1 but we recalibrate the luminosity evolution at higher redshifts to fol- 



l ow the extincted SF RD derived from the ultraviolet luminosity function of galaxies by 
(IBouwens et al.ll2009l ). This is simply the difference between their extinction corrected SFRD 
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derived using the UV slope and the unobscured SFR from direct measurements of the UV 
luminosity function. This differs from the earlier calculation at 2 > 1 in that it assumes 
that the IR luminosity density is flat out to 2 ~ 3 and then declines monotonically between 
3 < 2; < 6. Both these trends for the evolution of the IR LF are shown in Figure O We find 
that using the CEOl templates, we always exceed the FIR EBL at A > 300 /zm if we adopt 
a fiat IR luminosity density between 1 < 2; < 2.3 suggesting that the UV slope extinction 
estimates are too high. 



Another constraint at 2; > 1 comes from the radio stacking analysis of lCarilli et al.l (j2008[ ) 
who found that at 2 ~ 3, the ratio of the radio derived star-formation rate density to the 
UV-derived star-formation for the galaxies at the bright end of the UV luminosity function 
is 1.8±0.4. This implies that the dust obscured star-formation rate den sity is 0.8 times the 
star-formation rate derived from the observed UV values. Based on the iReddy et al.l (120061 ) 
measure of the UV luminosity function, this translates to a dust obscured star-formation 
rate of 0.02 yr^^ Mpc~'^. Due to the relatively shallow depth of the rest-frame UV 
observations in this field compared to the UDF, the fainter galaxies, a fraction of which are 
likely to be the heavily dust obscured LIRGs and ULIRGs, are missed. As a result, the 
derived value is most likely a lower limit. 

The range of possible constrai nts on the dust obs cured star-formation rate density de- 
rived above are consistent with the ICarilli et al.l (120081 ) measurement. We thereby adopt an 
= 1.3, ar, = —3 at 2 > 1 as the "best" model which is most consistent with the radio and 



the high redshift UV measurements of iBouwens et al.l (120091 ) . as well as the EBL constraints 
from FIRAS. 



3.2. 



Constraints on Dust Enshrouded Star- Formation at 2 > 1: Cooler FIR 

Templates 



Although the CEOl models have been shown to fit galaxy SEDs out to 2 ~ 1.3 (iMagnelli et al 



2009al ). there is evidence that due to the increased gas content of galaxies at higher redshifts, 
high redshift star-f orming galaxy SEID are different, mor e typical of scaled up low-luminosity 
starburst galaxies ( Pope et al. 2008, 2006 : Rigby et al.| [2008). Another possible reason for 
the evolving SED might be the decreasing metallicity which might be responsible for a change 
in grain size distribution. The discrepancy between the MIR to FIR ratios appears to be 
most striking at » 10^^ L©. Since these objects are only observed at 2 > 1.3, it is un- 
clear if this is a luminosity effect or a redshift effect or if there is a strong observational bias 
since we can evaluat e the properties of only the brightest, coolest galaxies at these redshifts 
(iMurphy et al.ll2009l ). We develop a revised library of SED in this paper which benefit from 
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the following observations: 



An extensive set of spectroscopic observations of in : 



LIRG and ULIRG regime with IRS flPope et al 



2008 



Tared luminous galaxies in the 



20091 . 12OO8I : iRigbv et aliboosl : lArmus et al.lboOTt lOesai et al-lboooh 



Murphy et al.ll2009l : ISiana et al. 



Deep 70 and 850 ^m imaging in the GOODS/FIDEL fields which co nstrains the long 



wavelength spectral energy distribution of the most luminous sources (jPope et al.ll2008 
Murphv et aPboOoh 



Stacking results at 70 /im as a function of 24 /im luminosity at z 
jMagneUi et al.ll2009al jbl) 



1 and z 



The data seem to suggest an increase in the MIR to FIR ratio of galaxies at high 
luminosities such that even the most luminous systems have large PAH equivalent width 
(Figure [7]). However, there are strong selection effects associated with detecting these objects 
at high redshift. As a result, it is not clear if the observed galaxies span the entire range of 
dust temperatures and associated far-infrared colors or are just the most extreme sources. 
We note that an extensi ve set of spect r oscop ic observations of a sample of bright 24 ^um source 
have been presented in iDasyra et al.l (120091 ). but these have not yet been incorporated into 
our database. We adopt this revised set of SED only at z > 1.3 since galaxies at z < 1.3 
appear to follow the CEOl templates quite well. 



In iMagnelli et al.l ( l2009bl ). sources in the redshift bin 1.3 < z < 2.3 were stacked at 
70 /im as a function of 24/im/(H-z) luminosity. It was shown that the CEOl templates 
overpredict the 70/im/(l+z) luminosity of the galaxies in the sense that the far- i n frared 
SE Ds of these o bjects are too warm. A similar trend was found in iMurphy et al.l (120091 ) 
and jPope et al.l ( l2008l ) although it was unclear if this is due to the fact that the sample was 
biased by the presence of submillimeter galaxies which might have cooler dust temperatures 
than average. We have modified the template SEDs to fit the observed trend between the 
24yum/(H-z) and 70 //m/(H-z) luminosity of the galaxies as discussed in Section 3. 

These new templates are shown in Figure [H] and compare favorably with the stacked 
spectral energy distribution of 10^^ Lq submillimet er galaxies at 2: ~ 2.5 for which we have 
high quality far-infrared constraints on the SED (jPope et al.ll2008l ). We emphasize that 
these SED are our first attempt to constrain the 2; ~ 1 — 2 far-infrared SED of galaxies and 
further work is needed to assess if these can be used as reliably indicators of the bolometric 
correction in 2; > 1.3 galaxies. 



We now evaluate constraints on the EBL using this set of cooler far-infrared templates 
that are in Figure [81 
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We use the CEOl templates for galaxies at z < 1.3 but for z > 1.3, we adopt these 
new cooler far-infrared templates. For z > 1, if we keep fixed at 1.3, then a/j must lie 
between —2.8 and —2.0 to be consistent with the EBL. Translating this to a dust obscured 
star-formation rate density using the canonical conversion factor of 1.71x10~^°Lq for 1 Mq 
we derive the upper limits on the dust obscured SFR density as shown in 
Figure [9l A comparison with Figure |5] indicates that these values are somewhat higher than 
the values derived using the CEOl templates but still overlap substantially. The similarity 
between these estimates and the estimates in the previous section suggest that the constraints 
on evolution at 2; > 1 are not extremely sensitive to slight changes in far-infrared SED of 
galaxies at high redshift. We thereby adopt = 1.3, an = —2.0 at z > 1 as our "maximal" 
dust obscured star-formation rate density model since it results in the highest star-formation 
rate density at 2; > 1 without violating the EBL values. 

A useful consiste ncy check on the m odel comes from the analysis of the 850 fim selected 
galaxy population by I Wall et al.l ( l2008l ). They found that submillimeter galaxies at 1.8 < 
z < 4 have a surface density of ~1800±600 deg~^ with the surface density at z < 1.8 being 
a factor of 5 larger. Based on the derived evolution at 2; > 1, we obtain a 2; > 1.8 surface 
density of 1100 deg~^ for the galaxy population with Sgso > 2mJy. Higher values for the 
surface density can be obtained simply by cutting off the faint end of the infrared luminosity 
function at ~1O^L0. Alth ough the "naaxima l" evolution derived here is ~2.5(j above the 
radio stacking constraint of ICarilli et al.l (120081 ). it indicates that the two models presented in 
this section and in Section 3.1 must straddle the entire range of 2; > 1 evolutionary histories 
for the IR luminosity function. These constraints on the dust obscured component of the 
star-formation rate density are compared with the unobscured star-formation rate density 
in Figure [TOl 

As we discussed in Section 3.1, we assess if it is possible to obtain consistency between 
these measurements and the UV slope estimates using the new templates. If the IR luminos- 
ity density is flat between 1 < z < 2.3 such that it is consisteri t wit h the lowest red s hift a t 
which a dust correct SFRD is measured by iReddy et al.l (120061 ) and iBouwens et al. J2009h . 
then the decline of the dust obscured star-formation rate density with increasing redshift 
must necessarily be steeper than measured from the UV to avoid violating the FIRAS mea- 
surements. We find that if the decline starts at 2; = 2.3, then for ai = 0, must lie in 
the range —6 < ao < —3.4. This factor of ~4 overestimate of the UV SFRD at 2; > 2.3 is 
most likely due to significant contribut ion of nebula r emis si on to the SEP of a young stellar 
population, an effect first noticed by iReddy et al.l (120061 ): ISiana et al.l (120081). They con- 



cluded that for galaxies with <10 Myr stellar populations, the UV slope and iMeurer et al. 
(I1999I ) relation overestimates the reddening partly due to the nebular emission and partly 
because the intrinsic colors of the stellar population are significantly different from the nom- 
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inal ~100 Myr old stellar population template adopted as a benchmark for the intrinsic SED 
of a galaxy. 

We next translate our constraints on the infrared luminosity density at ^ > 2 to an 
effective ultraviolet extinction for comparison with the extinction derived from the UV slope. 



3.3. Evolution of Dust Extinction 



It has previously been suggested by ICarilli et al.l (|2008[ ) th at the dust correct ions to 
the UV based star-formation rate density might be overestimates. ICarilli et al.l ( 120081 ) based 
their argument on the radio stacking of UV selected LBGs and found that at z ~ 3, the ratio 
of the radio derived star-formation rate density to the UV-derived star-formation is 1.8±0.4. 
If this ratio is attributed to dust extinction, then the extinction value is a factor of ~ 2 — 3 
lower than previous estimates which have found extinction corrections of ~3-5. 

We find that the range of values we derive for the dust obscur ed star-formation i n 
the previous subsection are in agreement with the values derived by ICarilli et al.l ( 120081 ). 
Thus, we conclude that the UV slope derived star-formation rate values are overestimates 
and the radio and IR are giving consistent results. Moreover this decline in SFRD between 



1 < z < 2.3 is in agreement with the estimates derived from 70 /im stacking in lMagnelli et al. 
(l2009bh . 



Based on these strong constraints from the EEL on the fraction of dust obscured star- 
formation at 2 > 1, we can utilize the measured UV luminosity density at these redshifts to 
predict the average extinction. We adopt the reddening-un corrected UV luminosity density 
down to 0.04L^ ._Trv^z=,ci measured bv iBouwens et al.l (120091 ) between 3.5 < z < 6 and those 



measured by iReddy et al.l (120061 ) at 1.9 < z < 3.4. We adopt the "best" constraint on 
the dust-obscured star- formation rate density derived using the modified CEOl templates 
in section 3.1 with ao at z > 1 of —3.0 although we provide uncertainties which span 
the range of possib le FIR values derive d from the analysis in Section 3.1 and 3.2. We use 
the prescription of iMeurer et al.l ( 119991 ) to convert the ratio of far-infrared to ultraviolet 
luminosity density to an ultraviolet extinction. This requires converting our total infrared 
(8 — 1000 um) luminos ity density to the far-infrared (40 — 120 /im) luminosity density used 
by lMeurer et al.l ( 1l999l ). As discussed in CEOl, this involves dividing the infrared luminosity 
by a factor of 1.92 for the galaxies in the IRAS Bright Galaxy Sample and so we apply the 
same ratio here. We have verified that this is a reasonable ratio to use even for the newer 
templates presented here by integrating the far-infrared SED over the relevant wavelength 
ranges although the conversion varies from 1.5 — 2.4 depending on the color temperature of 
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the far- infrared emission. 



Using the iMeurer et ahl ( 1999 ) conversion between the 1600 A extinction and the ratio 
of far-infrared luminosity to ultraviolet luminosity, we find that the upper limit to the dust 
obscuration at z > 1 imposed by the far-infrared EBL, results in extinction values that 
are systematically smaller than those derived from the UV slope. The discrepancy is the 
largest in the 2 < z < 4 range where the difference between our extinction and the UV slope 
extinction is a factor of ~2-3 (Figure [TT]) . There are two likely explanations for this. One 
possibility is that the increased fractional contribution due to nebular emission, which would 
have a relatively red SED, is making the UV slopes steeper. The other possibility is that 
the dust extinction law is changing due to the redshift evolution of dust grain sizes. For 
example, a Small Magellanic Cloud extinction law would decrease the ratio of the total to 
unobscured star formation rate estimate for a fixed value of the UV slope compared to the 
LMC or starburst extinction law which has been used to calibrate this relation (H. Shim, 
personal communication). 



Reddy fc Steidell (120091 ) have previously argued that at z ~ 2 — 3, the bulk of the 



star-formation occurs in faint, blue galaxies due to the steep faint end slope of the UV 
luminosity function. This is certainly plausible with the constraints on the dust obscured 
star-formation presented here. However, due to the uncertainty associated with converting 
UV-slope to extinction, errors are introduced in estimating dust obscured star-formation 
rate densities from the UV data. Current values derived from the UV are overestimates and 



should be reduced although the net effect on the comoving SFR measured by lReddy fc Steidel 
( I2OO9I ) is small. This would be in direct contrast to previous estimates of dust obscuration 
at 2; ~ 2, which indicate that ULIRGs dominate the comoving star- formation rate density. 
Furthermore, this is also in sharp contrast to the estimates at 2; < 1.3, where dust obscured 



infrared luminous galaxies a ccount for ~65% of the comoving SFR density (IMagnelli et al. 
2009al : ICharv fc ElbazI [200l[ ) . 



We also find good evidence for there to be an increase in average extinction with de- 
creasing redshift between 2 < z < 6. The extinction, which is correlated with dust content 
increases by a factor of 2 ove r this redshift range. Since the dust arises either from SNe or 
AGB stars ( iDwek et al.l 120071 ). the average dust content of the Universe must increase with 
stellar mass density. Over this redshift range, the stellar mass density has increased by a 
factor of ~ 4 suggesting that the dust content increases as the square root of the stellar mass 
density. In I Chary et al.l (j2007bl ). it was shown that the average metallicity in star-forming 
galaxies increases as stellar mass density to the power 0.69 ±0.1 7, over this redshift range. 
If metallicity were a proxy for dust content, the observed increased in extinction that we 
find in this paper, is consistent with the increase in metallicity that was found previously. 
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indicating that the two processes are unsurprisingly co-evolving. 



3.4. Comparison between IGL Estimates from Models vs Stacking 

We have used deep Spitzer imaging data at 24, 70 and 160 fim and submillimeter imaging 
from 250-850 /im to develop a reliable stacking procedure. In particular, we have been able 
to correct for the biases due to the clustering around bright detected sources in a model 
independent way. We have also used semi-empirical models for the evolution of the infrared 
luminosity function which is constrained by deep Spitzer imaging data at 24, 70 and 160 fim 
and the COBE measured EBL. These two independent lines of analysis should therefore given 
cross-consistent results on the estimates derived for the IGL. Figure [12] shows the values for 
the IGL derived from stacking compared with the EBL estimates derived from the model 
shown in Table 5. From the models, we have also estimated the fractional contribution to 
the EBL from LIRGs, ULIRGs and from bright 24 /im sources. We have used a different 
24 /im flux limit for each panel depending on the sensitivity of the data used for the stacking 
analysis. Contributions from LIRGs and ULIRGs have been derived by integrating the 
luminosity function between 10^^^^^ L© and > 10^^ Lq respectively. These values typically 
would depend on the scatter in the far-infrared spectral energy distribution of galaxies. 
For example, if a subsample of galaxies were substantially cooler than average, they would 
result in a higher fraction at longer wavelengths, but a lower fraction at shorter wavelengths. 
We have therefore calculated the range for both the modified CEOl templates described in 
Section 3.1 and the newer templates presented in Section 3.2. 

We find good agreement between the stacking values and the fractional contribution 
inferred from the evolution of the infrared luminosity function which gives us confidence 
that our conclusions are robust. There is seemingly a a discrepancy at 160 /xm since the 
stacked IGL is only ~20%, while stacking on all S24 > 55 /xJy sources should yield ~65%. 
This is due to the incompleteness from masking which as shown in Section 2.1.2 implies 
that we cannot stack below a 24 /im flux density of 275 /iJy. At a flux limit of 275 /iJy, the 
model reveals that we should obtain only about ~30% of the background, exactly what we 
derive from stacking. If we instead use the IGL value of 9.8nWm~^sr~^ that we get from 
simulating the 160 /im map with 24 /im sources down to 55/xJy, we find that we would have 
resolved ~70% of the EBL, in excellent agreement with the model. 

For future reference, we have also plotted the fractio n of the EBL that w ould be resolved 



by stacking on 850 /im and 1.2mm maps. We note that IWang et al.l (120061 ) have attempted 
stacking analysis on the 850 /im maps without subtracting the detections or masking. As 
discussed earlier, this results in severe overestimates of the stacking values and thereby the 
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derived IGL. A stacking analysis following similar techniques presented in this paper at these 
longer submillimeter wavelengths will be presented in Penner et al. (in preparation). 

We show in Figure [T3] the fractional contribution to the EBL at different wavelengths as 
a function of infrared luminosity. The largest contribution to the EBL arises from galaxies 
with luminosities between 1O^°~^^'^L0. Ultraluminous infrared galaxies with Lm > 10^^ 
contribute ~ 5 — 15% of the total EBL at any wavelength. In Figure [HI we plot the same 
quantity as a function of redshift which shows that the majority of the far-IR EBL comes 
from galaxies at z < 1.5. 

4. Predictions for Herschel 

Deep extragalactic surveys will be undertaken with the Herschel PACS and SPIRE 
instruments between 70 and 500 fxm as part of the Guaranteed Time surveys and the Guest 
Observer Legacy surveys. Using the derived constraints on the evolution of the far-infrared 
luminosity function of galaxies at z > 1, we can make predictions for the redshift distribution 
of far-infrared detected galaxies as a function of flux density and the fraction of the EBL 
resolved by these sources. 

4.1. Redshift Distribution 

We utilize both sets of templates, the modified CEOl and the new templates presented 
here to predict redshift distribution of sources as a function of far-infrared flux density. 
As discussed in Section 3.1 and 3.2, these correspond to our "best" model and "maximal" 
model. The redshift distributions are shown in Figure [151 The deepest Herschel surveys 
have the potential to detect dusty ultraluminous infrared galaxies out to 2 ~ 6 (solid curves), 
although the deep 500 iim survey with BLAST has potentially detected candidate objects in 
the 2; ~ 4 — 5 regime (dashed curves). Cross-identification and spectroscopic confirmation of 
these sources is likely to be difficult. Beyond ^ ~ 3, the PAH features get redshifted out of 
the 24 /im bandpass as a result of which the 24 /xm is a weak tracer of dusty galaxies at the 
highest redshifts. A deep 1.4 GHz radio survey of ~15 /xJy depth will be required to identify 
the counterparts of these galaxies, assuming the local radio-infrared correlation is valid at 
those redshifts. However, the surface density of 2: > 4 galaxies detected at 500 /um is likely 
to be low. We anticipate that there are ~2-9 such sources per GOODS field which covers 
about 0.046 deg2. 

Regardless, the shape of the redshift distribution of sources shows a peak around 2; ~ 1- 



- 28 - 



2 and then dec l ines r apidly. This is broadly similar to the redshift distribution presented in 
Valiante et al.l (|2009[ ) for Herschel surveys which is unsurprising since their star-formation 
rate density is quite similar to that shown in Figure 7 and 8. Our number counts appear 
to be quite similar t o their estimates as well. At a lOO^um flux density limit of 1.5 mJy, 
Valiante et al.l (120091 ) predict >10^ sources per square degree while we predict 1.4 x 10^ deg~^. 
At longer wave lengths though , signifi cant discrepancies creep in. Above a 350/im fiux density 
limit of 6 mJy, IValiante et al.l ( l2009l ) predict ~10'' deg"^ while we predict 2500 — 3500 deg~^ 
where the range of values arise from our "best" model and "maximal" model. In comparison, 
Patanchon et al. (2009) derive a source density of ~2900 deg~^ above 6 mJ y from a P(D) 



analys is of the BLAST maps. Similarly at a 500/im fiux density limit of 5 mJy. IValiante et al. 
J2009h predict 5000 deg-^ while we expect 1000 - 1500 deg '^ Patanchon obtai n ~2000 
deg~^ at the same limit. The origin of this difference between IValiante et al.l (120091 ) and our 



analysis is likely because their fits pass through the upper range of values in the submillimeter 
source counts since they adopted a fiat evolution for the infrared luminosity density between 
2; ~ 1 — 2 while our estimates favor a decline at 2; > 1. Their estimate should have resulted 
in an overproduction of the EBL although it is unclear from their Figure 16 whether their 
measurement is in agreement with the FIRAS limits. 

We conclude that both our estimates for the redshift distribution are in agreement but 
that we require a slightly steeper evolution for the z > 1 dust obscured star-formation rate 
to ensure consistency with the FIR EBL which results in lower source counts at the longest 
Herschel wavelengths. 



4.2. Fraction of EBL Resolved 

Direct detections with Herschel, even with the deepest extragalactic surveys such as 
GOODS-H, will resolve only a fraction of the total FIR EBL as shown in Figure [161 This 
fraction is 80%, 48%, 39%, 29% and 12% from sources brighter than fiux density limits 
of 0.6, 5, 4.2, 6 and 5 mJy at the following wavelengths: 100, 170, 250, 350 and 500 /im. 
However, the superior spatial resolution of Herschel compared to either Spitzer or BLAST 
will alleviate the effect of confusion noise and allow stacking to be done down to fainter limits 
than presented here. As shown in Figure [121 stacking exclusively on 24 /im galaxy samples 
is unlikely to result in significant improvement compared to the IGL estimates presented 
here since the fraction of the EBL that comes from the sample of 24 /im selected galaxies 
has been estimated to be ~ 40 — 60% and are already in agreement with our values derived 
from simulating the maps. Alternate possibilities such as stacking on near-infrared selected 
samples of sources will have to be investigated while ensuring that color cuts are used to 
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select the dusty sources from the sources that are not actively star-forming. We estimate 
that by stacking on 100s of dusty galaxies in differential flux bins identified through their 
optical to near-infrared colors, we can achieve noise values that are a factor of 10 below the 
direct survey limits quoted above. At these depths, the models indicate that 95%, 83%, 82%, 
71% and 68% of the background at 100, 170, 250, 350 and 500 /xm will be directly resolved 
which will be a significant improvement on the values quoted in Table |H 



5. Implications for Zodiacal Light Models and the Near-Infrared Extragalactic 

Background Light 



Given the consistency between our stacking analysis and the model f or the evolution 
of the infrared luminosity derived based on the work by iMagnelli et al.l (j2009al ). we feel 
confident to compare the EBL estimates at A < 140 /im derived from DIRBE data with the 
model estimates. As can be seen in Table 5, at these wavelengths, there is a significant 
discrepancy between the DIRBE estimates and the models. Constraints from the detection 
of TeV photons from blazars have excluded the possibility that the high value of the EBL at 
60 fim is due to an extragalactic source population. It is also clear that there is no unusual 
source population in the deep 70 iJ,m data that can account for the higher value of the 60 fim 
EBL. The most likely origin for the discrepancy is likely to be the zodiacal dust cloud. By 
comparing our model IGL intensities with these DIRBE EBL values, we attempt to constrain 
the origin of this excess emission. 

The interplanetary dust cloud, or zodiacal dust, is made up of dust from asteroid and 
comets. The thermal component of emission from the zodiacal cloud which dominates at 
A > 3.5 ^m and the sunlight scattered off the cloud at shorter wavelengths are the primary 
sources of foreground emission at infrared wavelengths. Models that have been developed 
for the zodiacal light have utilized the time varying component of the emissio n as the Earth 
revolves around the Sun ( iReach et al.l l2003l : iKelsall et al.l Il998l : IWrightl Il998l ) . These three 
dimensional models broadly indicate that the dust has a temperature of ~250K with a 
characteristic distance of ~1 Astronomical Unit (AU) from the Sun, although the exact values 
are a function of ecliptic latitude. Using the temporal variation of the emissi on however, 



provi des poor constraints on the isotropic component of the zodiacal emission ( iDwek et al. 
2005h . 



Wrightl ( 119971 ) tried to constrain any additional component of zodiacal light by applying 



the "very strong no zody" condition which states that the high Galactic latitude 25 //m 
intensity, where the zodiacal emission peaks, should be constant and zero in one of the 
DIRBE observations. Although this increased the contribution of zodiacal light emission to 
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the diffuse intensity at both near- and far-infrared wavelengths, the resultant extragalactic 
background light at near-infrared wavelengths was still substantia lly in excess of the intensit y 



obtained by adding the light from individually detected galaxies ( iLevenson fc Wrightll2008l ). 
For example, at 3.6 /im, corresponding to a minimum in the intensity of the z odiacal light 
emission, the integrated galaxy light is thought to be ~ 5.7 — 9nWm~^sr~^ (IFazio et al. 



2004J : ILevenson fc Wrightl 120081) while the EB L, using the "very strong no zody" condition 
is 13.3±2.8nWm~^ sr~^ (jLevenson et al. 2007 ). The discrepancy at shorter wavelengths i.e. 



1.2 and 2.2 /im, is larger, which may be a result of the increas ing contribution from the 



zodiacal dust scattering of direct sunlight ( ICambresy et al.ll200ll ). 



The discrepancy at near-infrared wavelengths that is discussed in the literature, com- 
bined with the difference between the EBL and the IGL at far-infrared wavelengths suggests 
a common origin for the emission. In Figure [T71 we plot the difference between the IGL 
and the EBL at both near- and far- infrared wavelengths. The discrepancy is maximum at 
60 /im, with a decreasing contribution out to 240 /im. Although, the significance of the dis- 
crepancy at any one wavelength varies from 2.9cr at 60 /im to 1.1a at 240 /im, the combined 
intensity between 60—240 /im is significant at the 3.2a" level after including the systematic 
error associated with the evolution of the far- infrared luminosity function. As a result, this 
discrepancy is considered significant enough to warrant further investigation. 

As shown in Figure [T71 the spectrum of this excess emission is inconsistent with the 
spectrum of the thermal emission in the zodiacal light models. The spectrum is substantially 
redder than the zodiacal emission between 60 — 240 /im. Even if we fit only the excess at 
the most significant data point, at 60 /im, by zodiacal emission, it falls short of explaining 
the difference between the IGL and EBL at both 1.2 and 2.2 /im although the discrepancy 
at each of the wavelengths is only slightly larger than la. The chi-square fit to the residual 
using the standard zodiacal model template results in a ~300 suggesting that incorrect 
zodiacal light subtraction cannot be responsible for the entirety of this residual emission. 

We estimate the temperature of the blackbody corresponding to this extra component 
of emission. We find that the excess emission can be well fit by a blackbody of 53 K with a 
95% confidence interval of 16 K. Due to the low significance of the discrepant emission, we 
do not attempt to constrain the emissivity of the blackbody although we note that the best 
fit with emissivity as a free parameter is consistent with an emissivity index of zero. The 
total thermal emission from this component is a mere 0.4% of the main zodiacal emission 
in even a dark field like the GOODS field and therefore could easily have been missed. The 
best fit temperature implies that this component is most likely in the outer solar system and 
probably arises from either the Kuiper Belt region or from outgassing of comets as they travel 
past their perihelion point along their orbit (Chary & Schlichting, in prep.). The existence 
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of such a component has previously been suggested by ISternI (119961 ) although the search for 
such a component in DIRBE data has previously been limited to low ec liptic latitudes where 
the component is thought to be most significant (IBackman et al.lll995l ). 



We now assess what the albedo of this component has to be to account for the dis- 
crepancy between the integrated galaxy light and the EBL at near-infrared wavelengths. 
At 1.2 /xm, we adopt the EBL values of ICambresv et al.l ( 120011 ). while at 2.2 and 3.6 /im, 



we adopt the EBL estimates of iGoriian et al. 



tood ). 



which are consistent with the more 



recent estimates of Levenson fc Wrightl (20081). For th e integrated galaxy light, at 1.2 and 



2.2 ^m, we adopt the values o f ICambresv et al.l (I20011). At 3 . 6 )um, we adopt the mean of 
the estimates by Fazio et al. (j2004 ) and Levenson &: Wright ( 2008 ). Table 5 summarizes 
the measurements. We find that if we assume the Kuiper Belt dust has the same grain size 
distribution as the zodiacal dust cloud, the albedo of the grains has to be ~4-5 times greater 
than the zodiacal light to reproduce the discrepancy in the NIR background, assuming that 
the grains are at ~ 40 AU. Since the grains in the zodiacal dust cloud are thought to have 
a NIR albedo of ~0.2 jKelsall et al.lll998[ IWrightl ^9981 ) . the high implied albedo of >0.8 
for this component of emission is consistent with grains that have an ice mantle. Remark- 
ably, with a change in a single parameter, the albedo, we are able to simultaneously fit the 
difference between the EBL and IGL at all three NIR wavelengths. 

We thereby conclude that there is evidence for the existence of a distribution of ice 
mantle grains at high ecliptic latitudes at a distance of ~40 AU, that is responsible for the 
discrepancy between the integrated light from galaxies and the cosmic infrared background 
light at both near- and far-infrared wavelengths. The total thermal emission from this 
component is >25nWm~^ sr~^, which is a small fraction of the zodiacal light emission even 
in high ecliptic latitude regions of sky. The relative contributions of the different components 
are shown in Figure [13 Future work will reveal the origin and physical properties of grains 
contributing to this component of emission. 



6. Conclusions 

We have utilized data from the deepest far-infrared surveys currently undertaken to eval- 
uate the fraction of the extragalactic background light that has been resolved into individual 
galaxies. By using a combination of direct detections and stacking on the imaging data, we 
demonstrate that more than ~50% of the background currently remains unresolved at all 
wavelengths between 100 — 500 /im. The 70 /im IGL, after large completeness corrections 
however comprises ~75% of the extragalactic background light derived from models for the 
evolution of the galaxy luminosity function. We also demonstrate that previous estimates 
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of the integrated galaxy light using stacking have significantly overestimated the values by 
multiple counting of the fiux from clustered sources and provide steps through which these 
can be corrected. 

We have then, extrapolated measurements of the infrared luminosity function measured 
at z < 1.2 to higher redshifts and use the resolved fraction of the background along with 
the DIRBE/FIRAS measured EBL to derive our best estimates of the dust obscured star- 
formation rate at z > 1. We find that the UV slope appears to overestimate the dust 
extinction at 2; > 2, most likely due to the significant contribution of nebular emission to the 
broadband ultraviolet photometry. We provide average measures of the extinction at these 
redshifts and find that the average extinction increases with decreasing redshift. If translated 
to an increase in dust content, this is consistent with the rate of increase in metallicity with 
redshift, suggesting that the two processes are co-evolving. We also quantify the fractional 
contribution to the far-infrared EBL as a function of both redshift and luminosity. We find 
that ultraluminous infrared galaxies contribute ~5% of the total EBL and that the bulk of 
the star-formation occurs in galaxies with luminosities ~10^°~^^'^ L0. Finally, we use the 
difference between the EBL derived from our models and the DIRBE measured values to 
provide tentative evidence for a contribution from ice mantle dust in the outer solar system. 
This dust appears to have a total intensity of ~25nWm~^sr~^ at high echptic latitudes and 
can account for the discrepancy between the IGL and EBL at both near- and far-infrared 
wavelengths. 
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7. APPENDIX 



We compare the flux es of sources measured in beam-convolved BLAST maps as mea- 
sured by iDye et al.l ( 120091 ) with the flux measured through aperture photometry in the un- 
convolved maps. We flnd that the aperture flux show a larger dispersion at fainter fluxes 
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where the fractional contribution of other sources becomes increasingly significant (Figure 
[T9l) . As a result, we choose to stack the pixel flux in the beam-convolved maps for the 
BLAST data. 

We also compare the distribution of BLAST pixel value before and after sources are 
subtracted in Figure [201 Subtracting the individually detected sources removes the positive 
skew in pixel values and results in a map which is dominated by the sum of confusion and 
instrumental noise. We note that the confusion noise in the submillimeter data is symmetric 
around zero since the beam has significant negative sidelobes while in the far-infrared data, 
confusion noise will only have a positive skew. 

We also assess how much of a role the clustering of galaxies affect the stacks by measuring 
the surface density of 24 /zm detected sources in the vicinity of a bright 500 /im detected 
source and comparing that to the average surface density of 24 detected sources over the 
entire field. Figure [21] illustrates that the bright submillimeter galaxies have a factor of two 
overdensity of 24 /im detected galaxies within a radius R compared to the average surface 
density of such sources. This suggests that the 500 fim sources may be detected because it 
is the combined flux of the individual sources. It also indicates that when stacking at the 
positions of 24 /xm sources, the 500 fim flux from the source will thus be counted twice unless 
the PSF is fitted and subtracted. Stacking on maps with sources still in it will result in a 
positive bias to the stacked results. 

REFERENCES 

Armus, L., Charmandaris, V., Bernard-Salas, J., Spoon, H. W. W., Marshall, J. A., Higdon, 
S. J. U., Desai, V., Teplitz, H. I., Hao, L., Devost, D., Brandl, B. R., Wu, Y., Sloan, 
G. C, Soifer, B. T., Houck, J. R., & Herter, T. L. 2007, ApJ, 656, 148 

Backman, D. E., Dasgupta, A., & Stencel, R. E. 1995, ApJ, 450, L35+ 

Bernstein, R. A. 2007, ApJ, 666, 663 

Bethermin, M., Dole, H., Beelen, A. & Aussel, H., 2010, A&A, in press (astro-ph/1001.0896) 
Bouwens, R. et al. 2009, ApJ 

Cambresy, L., Reach, W. T., Beichman, C. A., & Jarrett, T. H. 2001, ApJ, 555, 563 

Caputi, K. I., Lagache, G., Yan, L., Dole, H., Bavouzet, N., Le Floc'h, E., Choi, P. I., Helou, 
G., & Reddy, N. 2007, ApJ, 660, 97 



-34- 



Carilli, C. L., Lee, N., Capak, P., Schinnerer, E., Lee, K.-S., McCraken, H., Yun, M. S., 
Scoville, N., Smolcic, V., Giavalisco, M., Datta, A., Taniguchi, Y., & Urry, C. M. 

2008, ApJ, 689, 883 

Chary, R., 2007a, ASP Conference Proceedings, 380, 375 
Chary, R., Berger, E., & Cowie, L. 2007, ApJ, 671, 272 
Chary, R. & Elbaz, D. 2001, ApJ, 556, 562 
Dale, D. A. & Helou, G. 2002, ApJ, 576, 159 

Dasyra, K. M., Yan, L., Helou, G., Sajina, A., Fadda, D., Zamojski, M., Armus, L., Draine, 

B. , & Prayer, D. 2009, ArXiv e-prints 

Desai, V., Soifer, B. T., Dey, A., Le Floc'h, E., Armus, L., Brand, K., Brown, M. J. I., 
Brodwin, M., Jannuzi, B. T., Houck, J. R., Weedman, D. W., Ashby, M. L. N., 
Gonzalez, A., Huang, J., Smith, H. A., Teplitz, H., Willner, S. P., & Melbourne, J. 

2009, ArXiv e-prints 

Devlin, M. J., Ade, P. A. R., Aretxaga, I., Bock, J. J., Chapin, E. L., Griffin, M., Gundersen, 
J. O., Halpern, M., Hargrave, P. C, Hughes, D. H., Klein, J., Marsden, G., Martin, 
P. G., Mauskopf, P., Moncelsi, L., Netterfield, C. B., Ngo, H., Olmi, L., Pascale, E., 
Patanchon, G., Rex, M., Scott, D., Scmisch, C, Thomas, N., Truch, M. D. P., Tucker, 

C, Tucker, G. S., Viero, M. P., & Wiebe, D. V. 2009, Nature, 458, 737 

Dole, H., Lagache, G., Puget, J.-L., Caputi, K. I., Fernandez-Conde, N., Le Floc'h, E., 
Papovich, C, Perez- Gonzalez, P. G., Rieke, G. H., & Blaylock, M. 2006, A&A, 451, 
417 

Dwek, E., Arendt, R. G., & Krennrich, F. 2005, ApJ, 635, 784 
Dwek, E., GaUiano, F., & Jones, A. P. 2007, ApJ, 662, 927 

Dye, S., Ade, P. A. R., Bock, J. J., Chapin, E. L., Devlin, M. J., Dunlop, J. S., Eales, S. A., 
Griffin, M., Gundersen, J. O., Halpern, M., Hargrave, P. C, Hughes, D. H., Klein, 
J., Magnelh, B., Marsden, G., Mauskopf, P., Moncelsi, L., Netterfield, C. B., Olmi, 

L., Pascale, E., Patanchon, G., Rex, M., Scott, D., Semisch, C, Targett, T., Thomas, 
N., Truch, M. D. P., Tucker, C, Tucker, G. S., Viero, M. P., & Wiebe, D. V. 2009, 
ArXiv e-prints 

Elbaz, D., Cesarsky, C. J., Chanial, P., Aussel, H., Franceschini, A., Fadda, D., & Chary, 
R. R. 2002, A&A, 384, 848 



-35 - 



Engelbracht, C. W., Rieke, G. H., Gordon, K. D., Smith, J.-D. T., Werner, M. W., Mous- 
takas, J., Willmer, C. N. A., & Vanzi, L. 2008, ApJ, 678, 804 

Fazio, G. G., Ashby, M. L. N., Barmby, P., Hora, J. L., Huang, J.-S., Pahre, M. A., Wang, 
Z., Willner, S. P., Arendt, R. G., Moseley, S. H., Brodwin, M., Eisenhardt, P., Stern, 

D. , ToUestrup, E. V., & Wright, E. L. 2004, ApJS, 154, 39 

Finkbeiner, D. P., Davis, M., & Schlegel, D. J. 2000, ApJ, 544, 81 

Fixsen, D. J., Dwek, E., Mather, J. C, Bennett, C. L., & Shafer, R. A. 1998, ApJ, 508, 123 

Frayer, D. T., Huynh, M. T., Chary, R., Dickinson, M., Elbaz, D., Fadda, D., Surace, J. A., 
Tephtz, H. I., Yan, L., & Mobasher, B. 2006, ApJ, 647, L9 

Gilh, R., Daddi, E., Chary, R., Dickinson, M., Elbaz, D., Giavahsco, M., Kitzbichler, M., 
Stern, D., & VanzeUa, E. 2007, A&A, 475, 83 

Gorjian, V., Wright, E. L., & Chary, R. R. 2000, ApJ, 536, 550 

Hauser, M. C, Arendt, R. G., Kelsall, T., Dwek, E., Odegard, N., Weiland, J. L., Freuden- 
reich, H. T., Reach, W. T., Silverberg, R. F., Moseley, S. H., Pei, Y. C, Lubin, P., 
Mather, J. C, Shafer, R. A., Smoot, G. F., Weiss, R., Wilkinson, D. T., & Wright, 

E. L. 1998, ApJ, 508, 25 

Kelsall, T., Weiland, J. L., Franz, B. A., Reach, W. T., Arendt, R. C, Dwek, E., Freuden- 
reich, H. T., Hauser, M. G., Moseley, S. H., Odegard, N. P., Silverberg, R. F., & 
Wright, E. L. 1998, ApJ, 508, 44 

Lagache, G., Abergel, A., Boulanger, F., Desert, F. X., & Puget, J.-L. 1999, A&A, 344, 322 

Le Floc'h, E., Papovich, C, Dole, H., Bell, E. F., Lagache, G., Rieke, G. H., Egami, E., 
Perez-Gonzalez, P. G., Alonso-Herrero, A., Rieke, M. J., Blaylock, M., Engelbracht, 
C. W., Gordon, K. D., Hines, D. C, Misselt, K. A., Morrison, J. E., & Mould, J. 
2005, ApJ, 632, 169 

Levenson, L. R. & Wright, E. L. 2008, ApJ, 683, 585 

Levenson, L. R., Wright, E. L., & Johnson, B. D. 2007, ApJ, 666, 34 

MagneUi, B., Elbaz, D., Chary, R. R., Dickinson, M., Le Borgne, D., Frayer, D. T., & 
Willmer, C. N. A. 2009a, A&A, 496, 57 

— . 2009b, A&A 



-36- 



Marsden, G., Ade, P. A. R., Bock, J. J., Chapin, E. L., Devlin, M. J., Dicker, S. R., 
Griffin, M., Gundersen, J. O., Halpern, M., Hargrave, P. C., Hughes, D. H., Klein, 
J., Mauskopf, P., Magnelli, B., Moncelsi, L., Netterfield, C. B., Ngo, H., Olmi, L., 
Pascale, E., Patanchon, G., Rex, M., Scott, D., Semisch, C., Thomas, N., Truch, 
M. D. P., Tucker, C., Tucker, G. S., Viero, M. P., & Wiebe, D. V. 2009, ArXiv 
e-prints 

Meurer, G. R., Heckman, T. M., & Calzetti, D. 1999, ApJ, 521, 64 

Murphy, E. J., Ghary, R.-R., Alexander, D. M., Dickinson, M., Magnelli, B., Morrison, G., 
Pope, A., & Teplitz, H. 1. 2009, ApJ, 698, 1380 

Pattanchon, G., et al. 2009, ApJ, 707, 1750 

Pascale, E., Ade, P. A. R., Bock, J. J., Chapin, E. L., Devlin, M. J., Dye, S., Eales, S. A., 
Griffin, M., Gundersen, J. O., Halpern, M., Hargrave, P. C., Hughes, D. H., Klein, J., 
Marsden, G., Mauskopf, P., Moncelsi, L., Netterfield, C. B., Olmi, L., Patanchon, G., 
Rex, M., Scott, D., Semisch, C., Thomas, N., Truch, M. D. P., Tucker, G., Tucker, 
G. S., Viero, M. P., & Wiebe, D. V. 2009, ArXiv e-prints 

Pope, A., Chary, R.-R., Alexander, D. M., Armus, L., Dickinson, M., Elbaz, D., Prayer, D., 
Scott, D., & Teplitz, H. 2008, ApJ, 675, 1171 

Pope, A., Scott, D., Dickinson, M., Chary, R.-R., Morrison, G., Borys, C., Sajina, A., 
Alexander, D. M., Daddi, E., Prayer, D., MacDonald, E., & Stern, D. 2006, MNRAS, 
370, 1185 

Puget, J.-L., Abergel, A., Bernard, J. -P., Boulanger, P., Burton, W. B., Desert, F.-X., & 
Hartmann, D. 1996, A&A, 308, L5+ 

Reach, W. T., Morris, P., Boulanger, F., & Okumura, K. 2003, Icarus, 164, 384 

Reddy, N. A. & Steidel, C. C. 2009, ApJ, 692, 778 

Reddy, N. A., Steidel, C. C., Fadda, D., Yan, L., Pettini, M., Shapley, A. E., Erb, D. K., & 
Adelberger, K. L. 2006, ApJ, 644, 792 

Rieke, G. H., Alonso-Herrero, A., Weiner, B. J., Perez-Gonzalez, P. G., Blaylock, M., Donley, 
J. L., & Marcillac, D. 2009, ApJ, 692, 556 

Rieke, G. H., Young, E. T., Engelbracht, C. W., Kelly, D. M., Low, F. J., Haller, E. E., 
Beeman, J. W., Gordon, K. D., Stansberry, J. A., Misselt, K. A., Cadien, J., Morrison, 
J. E., Rivhs, G., Latter, W. B., Noriega-Crespo, A., Padgett, D. L., Stapelfeldt, K. R., 



-37- 



Hines, D. C, Egami, E., MuzeroUe, J., Alonso-Herrero, A., Blaylock, M., Dole, H., 
Hinz, J. L., Le Floc'h, E., Papovich, C, Pcrez-Gonzalcz, P. G., Smith, P. S., Su, 
K. Y. L., Bennett, L., Prayer, D. T., Henderson, D., Lu, N., Masci, P., Pesenson, 
M., Rebull, L., Rho, J., Keene, J., Stolovy, S., Wachter, S., Wheaton, W., Werner, 
M. W., & Richards, P. L. 2004, ApJS, 154, 25 

Rigby, J. R., Marcillac, D., Egami, E., Rieke, G. H., Richard, J., Kneib, J.-R, Fadda, D., 
Willmer, C. N. A., Borys, C., van der Werf, P. P., Perez- Gonzalez, P. G., Knudsen, 
K. K., & Papovich, C. 2008, ApJ, 675, 262 

Sanders, D. B., Mazzarella, J. M., Kim, D.-C., Surace, J. A., & Soifer, B. T. 2003, AJ, 126, 
1607 

Schiminovich, D., et al. 2005, ApJ, 619, L47 

Siana, B., Smail, I., Swinbank, A. M., Richard, J., Teplitz, H. I., Coppin, K. E. K., Ellis, 
R. S., Stark, D. P., Kneib, J.-R, & Edge, A. C. 2009, ApJ, 698, 1273 

Siana, B., Teplitz, H. I., Chary, R.-R., Colbert, J., & Prayer, D. T. 2008, ApJ, 689, 59 

Stern, S. A. 1996, A&A, 310, 999 

Treister, E. & Urry, C. M. 2006, ApJ, 652, L79 

Valiante, E., Lutz, D., Sturm, E., Genzel, R., & Chapin, E. 2009, ArXiv e-prints 
Wall, J. v.. Pope, A., & Scott, D. 2008, MNRAS, 383, 435 
Wang, W.-H., Cowie, L. L., & Barger, A. J. 2006, ApJ, 647, 74 

Wright, E. L. 1997, in Bulletin of the American Astronomical Society, Vol. 29, Bulletin of 
the American Astronomical Society, 1354 — h 

Wright, E. L. 1998, ApJ, 496, 1 



This preprint was prepared with the A AS macros v5.2. 



-38- 



Table 1: Description of the imaging data used for the stacking analysis. 



A (/um) 


PSF FWHM (") 


Pixel size (") 


R (pixels)'' 


Field 


S24//"Jy'' 


70 


18 


4 


4 


GOODS-N 


25 


70 








ECDF-S 


40 


160 


38 


8 


6.75 


EGS 


55 


250 


36 


10 


4.5 


ECDF-S 


40 


350 


42 


10 


5.5 


ECDF-S 


40 


500 


60 


10 


7.5 


ECDF-S 


40 



"R is the radius used for masking in the cutouts. 
''Flux hmit of 24 fim catalog used for stacking analysis. 
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Table 2: Contribution of detected far- infrared/sub millimeter sources to the IGL. 



Detection level 






uli, [nWm 


2 sr-1] 






70 yum 


160 /im 


250 /im 


350 /im 


500 /im 


6(7 


1.9 ±0.2 


l.liO.l 


0.40 ±0.01 


0.232 ±0.004 


0.103 ±0.002 


5(7 




1.5 ±0.2 


0.60 ±0.01 


0.356 ±0.006 


0.155 ±0.003 


4(7 




2.0 ±0.3 


0.97 ±0.01 


0.538 ± 0.008 


0.199 ±0.003 



Table 3. Contribution of 24 yum sources to the IGL as estimated from the data. 



Method 


70 fj.in 


160 fim 


ul^ [nW m^ sr 
250 Aim 


-i]b 

350 /j.m 


500 //m 


Dole et al. 2006 (> 60/iJy) 
Devlin et al. 2009 


5.9±0.9 


10.7±1.6 


7.9 ±0.5 


4.5 ± 0.3 


2.0 ±0.2 


No subtraction/no masking* 
Subtract 6a /no masking 
Subtract 5a /no masking 
Subtract Aa/no masking 


7.1±1.4 
5.9±1.1 (3.7) 


15.0±2.0 
13.2±1.7 (11.8) 


7.30 ± 0.46 
6.70 ± 0.40 (6.30) 
6.37 ±0.36 (5.77) 
5.84 ±0.32 (4.87) 


4.02 ± 0.24 
3.46 ±0.22 (3.23) 
3.28 ±0.20 (2.92) 
3.08 ±0.18 (2.54) 


1.86 ±0.14 
1.34 ±0.11 (1.24) 
1.26 ±0.11 (1.10) 
1.20 ±0.10 (1.00) 


Subtract Acr /masking 
Subtract Stj/masking 
Subtract 6(T/masking 
Subtract 6<T/masking'= 
No subtraction/masking 


4.2±0.7 (1.96) 
7.0±1.0 (4.8) 


3.3±0.4 (1.9) 


1.44 ± 0.09 (0.47) 
1.12 ±0.10 (0.52) 
0.97 ±0.11 (0.57) 

0.71 ±0.12 


0.72 ±0.04 (0.18) 
0.55 ±0.05 (0.19) 
0.45 ±0.05 (0.22) 

0.28 ±0.06 


0.223 ±0.016 (0.024) 
0.183 ±0.019 (0.028) 
0.138 ±0.021 (0.035) 

0.056 ± 0.029 



"These are derived by following exactly the same procedure as Dole et al. and Devlin et al. to enable a direct comparison with 
their values. 

''In cases where we have subtracted the detections out of the map we add their contribution after stacking. Values in brackets are 

from stacking only. 

"^Includes a completeness correction for faint galaxies which could not be stacked on. The completeness correction factor cannot be 
accurately determined for wavelengths other than 70 /xm. 
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Table 4: Estimate of the IGL through stacking simulations. The bias factor is the ratio of 
stacked flux with no masking to input flux as determined from simulations of the residual im- 
ages. Uncertainties on the bias factor reflect the range of templates used to convert between 
24 ^m and far-infrared flux and the noise in the maps. IGLcorr is the estimated contribu- 
tion to the background from 24 /im selected galaxies based on the simulations and has the 
contribution of individually detected sources added. Errors on IGLcorr include measurement 
errors as well as uncertainty in the calibration. 



Wavelength 


Bias factor 


IGLcorr 


fraction 






nW m^^ sr"'^ 




70 /im 


1.2 ±0.1 


5.3 ± 1.1 


~56% 


160 /im 


1.4 ±0.1 


9.8 ±2.0 


~ 70% 


250 /im 


1.16 ±0.05^ 


5.2 ±0.5 


~54% 


350 fxm 


1.14 ±0.05 


2.8 ±0.3 


~ 52% 


500 fim 


1.2 ±0.06 


1.0 ±0.2 


~43% 



The higher bias factor for 250 /im compared to 350 microns can be explained by the extended wings of the PSF at 250 /im. 



Table 5. Estimates of the EBL and IGL 



WcLvelengtn 


E13Ij Intensity 


rtererence 


ICjL Intensity 


rteterence 






11 V V 111 ol 




11 V V 111 oi 






1.2 


54 ± 16.8 


Cambresv et al. (2001) 


9 ± 2 


Cambresv et al. i 


[2001) 


2.2 


22.4 ± 6.0 


Gorlian et al. (2000) 


7.2 ± 1.6 


Cambresv et al. i 


[2001) 


3.6 


11.0 ± 3.3 


Gorlian et al. (2000) 


5.7-9.0 


Fazio et al. (2004): Levenson & Wrieht i 


[2008) 


60 


28.1 ± 7 


Finkbeiner et al. (2000) 


7.6 ± 0.25 


Model; This 


paper 


70 






9.5 ± 0.25 


Model; This 


paper 


100 


24.6 ± 8.4 


Finkbeiner et al. (2000) 


13.2 ± 0.5 


Model; This 


paper 


140, 160 


25 ± 7.0 


Hauser et al. (1998) 


14.1 ± 0.8 


Model; This 


paper 




12.6±4 


Fixsen et al. (1998) 










15.3±6.4 


Laeache et al. (1999) 








240 


14 ± 3.0 


Hauser et al. (1998) 


9.6 ± 0.9 


Model; This 


paper 




10.9±3 


Fixsen et al. (1998) 










11.4±1.9 


Lagache et al. (1999) 








350 


5.6 ± 1.7 


Fixsen et al. ( 1998) 


5.4 ± 0.8 


Model; This 


paper 


500 


2.4 ± 0.7 


Fixsen et al. (1998) 


2.3 ± 0.5 


Model; This 


paper 


850 


0.5 ± 0.15 


Fixsen et al. ( 1998) 


0.50 ± 0.14 


Model; This 


paper 


1200 


0.17 ± 0.05 


Fixsen et al. ( 1998) 


0.16 ± 0.05 


Model; This 


paper 
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Fig. 1. — Top Panel: The integrated galaxy light at 70 (left) and 160 /xm (right) from sources 
brighter than a particular 24 /im flux threshold, as estimated from Spitzer decta, using different 
stacking methodologies. The solid black is the best estimate from the current data. The 
EBL from DIRBE at 60 and 140 fj,m is shown as the sohd square, the EBL from FIRAS at 
160 /im as the solid star while the model estimates presented in this paper are shown as the 
empty square. When bright sources that are individually detected are not subtracted and 
the stacking is done without masking out the pixels that are included in the stack, the total 
IGL is significantly overestimated due to overcounting of the fiux (dotted line). Even after 
the sources are subtracted, stacking without masking results in the stacked intensity being 
biased high. After masking, large and uncertain completeness corrections are required to 
measure the IGL of sources below a 24 yum flux density value of 74 n^y and 275 /xJy in the 70 
and 160 fj,m images respectively. This is shown in the lower panels as the fraction of pixels 
attributed to sources in differential 24 fj,m flux bins, that are unmasked and can contribute 
to the stacked flux. Higher spatial resolution far-infrared data with Herschel will reduce the 
incompleteness due to masking and improve stacked estimates of the IGL from faint 24 fim 
sources. The solid black line in the left panel shows the 70 /im IGL after these completeness 
corrections arc applied. At 160 /i,m, the completeness correction can only be applied above 
275 /uJy and results in values similar to the solid line. As a result, we omit it for clarity. 
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Fig. 2. — Subtracting the submillimeter detections from the maps prior to stacking. Images 
are 27' across. Circles have a radius of 45" and indicate the positions of the 5a detections. 
The left panel shows the raw (decon) 250 /im map and the right panel is the same map with 
the detections removed. 
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Fig. 3. — Distribution of 24 iim flux densities of sources that can be stacked on at 250 //m 
after masking. The grey sohd curve shows ah secure 24 detections in the ECDF-S and 
the black dashed curve is the flux distribution for the sources we are able to stack in the 250 
/xm map with masking. The vertical dotted line shows the flux at which we are no longer 
complete in our stacking with masking. With only a few tens of sources we do not obtain a 
detection when stacking differential bins below = 600 /xJy. The dot-dash hne is the flux 
distribution of sources from stacking with masking on a map with same spatial resolution 
as 250 /im with Herschel; with 100s of sources per bin Herschel will be able to constrain the 
contribution to the IGL from sources with 5*24 = 20-600 /uJy. 



-46 - 



12 - 

Contribution from detections 
2Q FIRAS Stacked contribution no masking 

Stacked contribution with masking 




Nu sub 6a 5c 4.5c 4a 

NO MASKING 



WITH MASKING 



Contribution from detections 
PIRAS Stacked contribution no masking 

Stacked contribution with masking 



B 4 




No sub 6a 5c 4.5c 4a 

NO MASKING 



4o 4.5o 5c 6a No sub 

WITH MASKING 



3.0 
2.5 
2.0 



: FIRAS 



Contribution from detections 
Stacked contribution no masking 
Stacked contribution with masking 




Nu sub 6a 5c 4.5c 4a 

NO MASKING 



4o 4.5o 5c 6a No sub 

WITH MASKING 



Fig. 4. — Contribution of 24 /xm galaxies to the 250 (top left), 350 (top right) and 500 /im 
(bottom) submillimeter background. Light bars are the contribution of bright, individually 
detected submillimeter sources while dark bars are the contribution from stacking on 24 /im 
sources that are not associated with the bright submillimeter sources. The bottom axis 
indicated the signal-to-noise threshold down to which detected sources were subtracted from 
the maps. The bars on the left are for stacking without masking - these measurements will 
be biased high (as shown by the lower limit). The bars on the right are for stacking with 
masking - these measurements will be biased low (as shown by the lower limit). The blue 
data point is our best estimate for the background from all 24 /zm sources (S'24 > 40 /iJy) 
derived by simulating the stacking experiment. 
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Fig. 5. — Constraints on the z > 1 far-infrared luminosity density, and thereby dust obscured 
star-formation rate from the extragalactic background hght using the CEOl templates. The 
lu minosity density at z < 1 is constrained directly by far-infrared observations presented 



m 



Magnelli et al.l (l2009al ). The left panel shows different scenarios for the evolution of the 



far-infrared luminosity density while the right panel shows the resultant EBL (green lines) in 
comparison with observations. DIRBE estimates of the EBL are shown as blue stars while 
the red region is the FIR AS estimate of the EBL. Dashed lines are the early estimates of 
the far-infrared EBL from lPuget et al.l (119961 ). We allow for pure density evolution at 2; > 1, 
although density and luminosity evolution are degenerate as demonstrated in CEOl. The 
solid black star in the left panel is the radio stacking estimate on 2; ~ 3 LBGs by lCarilli et al. 



(120081). The do t -dash lines and points are estimates based on the UV-slope of galaxies from 



Bouwens et al.l (120091 ). We find that the dust obscured SFRD cannot be as hi gh as estimated 



by using the locally der ived reddening prescription on UV slope estimates ( jBouwens et al. 



2OO9I : iReddv et al.ll2006l ). 
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Fig. 6. — The left panel shows three possible estimates of the dust obscured star-formation 
rate density between < 2 < 6 (See text for det ails). The dot-dash lin es and points are 
estimates based on the UV-slope of galaxies from iBouwens et al.l (|2009[ ). The green lines 
in the right panel shows the inferred EBL if the SFR density were to evolve as in the left 
panel using the CEOl templates. Linestyles are preserved in the left and right panels in that 
the solid green line in the right panels corresponds to the solid black line in the left panel 
etc. The dust obscured SFRD estimated from the UV slopes appears to over predict the 
EBL by ~1.5cr at A > 350/xm suggesting that UV slopes overestimate dust obscuration at 
z > 2.5. The solid circles show the fraction of the FIR IGL which will be directly detected 
with Herschel; we will need to rely on stacking to resolve the rest of the background. 
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Fig. 7. — Ratio of the peak 7.7 fim PAH flux to the 9—11 /iin flux as a function of infrared 
luminosity. The ordinate is effectively a measure of the 7.7/im PAH equivalent width. Due 
to the difficulty in measuring the absolute continuum level in high redshift galaxies, we plot 
this ratio ins t ead o f an equivalent width. The solid black circles are local ULIRGs from 



Armus et al.l (120071). empty cyan diam onds are 0.5 < z < 2.5 infrared l umino us galaxies 



( iMurphy et al 



2009; 



Pope et al.l 120081) ■ empty stars are from ISiana et al.l (l2009[ l and aster- 



isks are from iRigby et al.l (l2008[ l. The solid black line shows the CEOl templates which 
underestimate PAH equivalent width among the most luminous sourc es due to the i n creas- 
ing AGN contribution in local ULIRGs. T he green dashed lin e are the iDale fc Heloul (|2002[ ) 
templates while the red dashed line are the iRieke et al.l (|2009[ ) SED models. Neither of these 
models reproduce the observed trend as a result of which we use have developed revised 
SEDs that follow the observational trend plotted as the black diamonds. 
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Fig. 8. — Comparison between the SEDs developed here (sohd hne) to fit the mid- infrared 
properties of z > 1.3 galaxies versus the CEOl templates (dotted line). The high redshift 
galaxies that have been observed appear to show higher PAH equivalent width and cooler 
dust temperature which motivated the developement of the new SEDs. However, it is unclear 
if this is a selection effect associated with the sensitivity of the existing far-infrared surveys. 
The IR luminosities of the plotted templates are 3.7x10^°, 2.8xl0^\ 2.8x10^^ ^nd lO^^L 



0- 



We note that the ne w 10 Lf^^ tenaplate is consistent with the stacked SED of submillimeter 
galaxies presented in iPope et al.l (120081 ) . 
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Fig. 9. — Same as Figure [5] but using the new cooler set of galaxy templates developed in 
this paper for objects at z > 1.3. The decline is marginally slower than for the modified 
CEOl templates. 
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Fig. 10. — The unobscured (blue hatched) and dust obscured (red hatched) star-formation 
rate density with cosmic time. The blue dotted hne is a s i mple polynomial fit to the ultravi 



olet measured SFR densities in ( ISchiminovich et al.ll2005t iReddy et al.ll2006l : iBouwens et al. 



20091 ) while the region around it spans the range of measured uncertainties. The red region 



is the obscured SFR density derived from our analysis. The co nstraints at z > 1 on the ob- 
scured SFR density are ambiguous. If the radio stacking point ( ICarilli et al.ll2008l ) is robust, 
then the lower limit of the region would be the best estimate of the obscured SFRD while 
the UV slope estimates would suggest the upper line. 
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Fig. 11. — Our best estimate of the I6OOA ultraviolet extinction (solid squares) derived from 
the models presented here along with the range of values (grey regions) that a re consistent 



with t he E BL. Empty circl e s show the extinction derived from the UV slope by iReddy et al. 



(I2OO6I ) and iBouwens et al.l (|2009( ) which would violate the FIRAS measured EBL estimates. 
The dust obscuration does increase by a factor of ~2 with decreasing redshift between 
2 < z < 6 which is not surprising since the dust content of the Universe must be correlated 
with the build up of the comoving stellar mass density which has increased by a factor of 
~4 between 2.3 < z < 5. 
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Fig. 12. — Fractional contribution to the EBL as a function of redshift with each panel 
representing a different wavelength. In each panel, we show the cumulative contribution of 
LIRGs, ULIRGs and galaxies above the 24 //m flux limit used to perform the stacking at 
each wavelength. For example, the GOODS-N data where the 70 fim stacking was done are 
sensitive down to ~25 /iJy while the ECDF-S data where the BLAST stacking was done is 
sensitive down to ~40 /xJy. The exact contribution of each population is dependent on the 
choice of FIR SED and so we show the plausible range for each quantity as the shaded region 
between the dotted and sohd lines. Also plotted are our estimates derived using stacking. 
The lower limit of the stacking point is the value derived by adding the individual detections 
to the stacked flux from other 24 fim sources, including masking. The upper limit is from 
adding the detection to the stacked flux using no masking. These limits are not measures of 
uncertainties, which are negligibly small but are provided to illustrate the difference in values 
if stacking is not done carefully. The best estimate IGL, including individual detections, 
stacking and with corrections for incompleteness and overcounting is plotted as the solid 

svmbnl Exr.pnt fnr Ififly/.Tn wHptp wp rarmnt stark bfilnw a. 9Aii,m fliTif nf 57.5 ;/,.Tv Hiip, tn 
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Fig. 13. — Fractional contribution to the EBL at different wavelengths as a function of IR 
luminosity, integrated over all redshifts. The dashed lines are for the CEOl templates as 
discussed in Section 3.1 while the solid lines are for the new templates presented in Figure 
6. LIRGs and ULIRGs (Lir > 10" Lq) contribute about ~35% of the total EBL at ~70 /im 
with the fraction increasing with increasing wavelength. 
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Redshift 

Fig. 14. — Fractional contribution to the EBL at different wavelengths as a function of 
redshift, integrated over all luminosities. The dashed lines are for the CEOl templates while 
the solid lines are for the new templates presented in Figure 6. More than ~70% of the EBL 
at A < 500 yum arises from galaxies at 2; < 1.5. 
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Fig. 15. — Predicted redshift distribution of deep far-infrared surveys between 100 and 
500 /im. The sohd hues are for a fainter flux density hmit (e.g. deep Herschel surveys such 
as GOODS-H) while the dashed lines are for a brighter flux density (e.g. current BLAST 
limits) at the same wavelength. The black lines show the estimates based on the "best" 
estimate of the IRLF at z > 1 as discussed in Section 3.1. The blue hues show the upper 
limits based on the "maximal" evolution in Section 3.2 which are consistent with the FIR 
EBL. 
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Fig. 16. — Fractional contribution to the EBL from sources above a specific flux density at 
different far-infrared wavelengths. The solid line is for the "best" model while the dashed 
line is for the "maximal" model. The deepest surveys with Herschel such as GOODS-H (P.I.: 
David Elbaz) will reach flux density limits of ~l-5 mJy. In order to resolve the bulk of the 
far-infrared/submillimeter EBL, we will have to stack hundreds of sources to achieve stacked 
flux density limits more than 10 times fainter than direct detection limits. 
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Fig. 17. — Left panel: The symbols show the discrepancy between the integrated galaxy light 
and the estimated value of the extragalactic background light at both near-infrared and far- 
infrared wavelengths. The solid hne shows a single component of ice mantle dust, with an 
albedo of ~1, at a distance of ~20-80 AU that can reproduce the observed discrepancy. The 
dashed line shows the intensity of the zodiacal dust cloud in GOODS-N (echptic latitude 
~ 60deg) while the dotted line is the same zodiacal light spectrum scaled down by 0.1. The 
observed discrepancy between the IGL and EBL is inconsistent with arising from the main 
cloud and requires a second component of emission from the outer solar system. Right panel: 
Chi-square contours from fitting the far-infrared emission by a blackbody. The different 
colors show the 90, 95 and 99% confidence levels while the best fit shown as the black square 
is /3 = 0, T = 53 K. 
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Fig. 18. — Estimates of the relative contributions of the zodiacal dust and integrated galaxy 
light in a typical extragalactic field. The latter would be isotropic and independent of choice 
of field. A tentative estimate of the contribution from isotropically distributed dust in the 
outer solar system is also labeled as cometary dust. Also shown is the intensity of the cosmic 
microwave background. 
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Fig. 19. — Aperture correction for photometry on 250 fim image cutouts. Taking the > 4cr 
sources in the central deep region of the BLAST ECDFS map, we compare the flux from 
the catalog to that measured from fitting the PSF to 9 x 9 cutouts and derive a negligible 
aperture correction of 1.02. At faint flux densities, the confusion from other sources increases 
resulting in a large scatter. 
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Fig. 20. — Pixel distribution for raw (light gray) and residual (black) 250 /im signal maps. 



-63- 



A 3 



V 



(N 







1 1 1 1 1 1 1 1 

- 
- 


1 1 1 1 1 


1 1 1 1 1 1 IV 1 1 1 1 1 1 1 1 

- 


o 


o 




o 




o 


o 




- 


___________ 




o o o = 


; O O 




O O O ; 


r O 




O \ 








III 


1 1 1 1 1 


/\ 1 1 1 1 


20 


30 


40 50 



S5oo(mJy) 



Fig. 21. — Ratio between the number of 24 //m sources within a radius R (see Table 1) of a 
BLAST 500 yum source and the average number of 24 fim sources expected within the same 
radius assuming a completely random distribution. The average is 2 (dashed line) and is 
independent of 500 /im flux density. This implies that when stacking at the positions of 
24 iim sources, the flux from individual 500 /xm sources will be double counted. 



